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Preface 


Despite  recent  advances  in  computer  hardware,  the  innovation  of  new,  efficient 
numerical  algorithms  still  yields  the  maximum  gain  in  analyzing  problems  once  thought 
to  be  too  large,  too  difficult,  or  too  expensive  to  solve.  I  believe  the  exponential 
characteristic  discrete  ordinates  scheme,  developed  here  in  slab  geometry,  demonstrates 
this  theme  extremely  well.  Expansion  of  this  scheme  to  more  complex  geometries  should 
produce  even  more  remarkable  results. 

I  would  like  to  thank  my  advisor,  LCDR  K.  Mathews,  Ph.D.  for  suggesting  this 
research  topic  and  providing  me  with  background  materials  and  existing  computer  codes. 
His  guidance  and  instinctive  numerical  methods  expertise  were  essential  in  completing 
this  project.  Also,  my  heartfelt  gratitude  goes  to  my  loving  wife  Patti,  who  faithfully 
supported  me,  and  to  God  for  granting  me  the  wisdom  and  patience  to  finish  this  project. 
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Abstract 

A  new  discrete  ordinates  spatial  quadrature  scheme  is  presented  for  solving  neutral 
particle  transport  problems.  This  new  scheme,  called  the  exponential  characteristic 
method,  is  developed  here  in  slab  geometry  with  isotropic  scattering.  This  method  uses  a 
characteristic  integration  of  the  Boltzmann  transport  equation  with  an  exponential 
function  as  the  assumed  form  of  the  source  distribution,  continuous  across  each  spatial 
cell.  The  exponential  source  function  is  constructed  to  globally  conserve  zeroth  and  first 
spatial  source  moments  and  is  non-negative.  Characteristic  integration  ensures 
non-negative  fluxes  and  flux  moments.  Numerical  testing  indicates  that  convergence  of 
the  exponential  characteristic  scheme  is  fourth  order  in  the  limit  of  vanishingly  thin  cells. 

Highly  accurate  solutions  to  optically  thick  problems  can  result  using  this  scheme 
with  veiy  coarse  meshes.  Comparing  accuracy  and  computational  cost  with  existing 
spatial  quadrature  schemes  (diamond  difference,  linear  discontinuous,  linear 
characteristic,  linear  adaptive,  etc.),  the  exponential  characteristic  scheme  typically 
performed  best.  This  scheme  is  expected  to  be  expandable  to  two  dimensions  in  a 
straight  forward  manner.  Due  to  the  high  accuracies  achievable  using  coarse  meshes,  this 
scheme  may  allow  researchers  to  obtain  solutions  to  transport  problems  once  thought  too 
large  or  too  difficult  to  be  adequately  solved  on  conventional  computer  systems. 


I. 


Neutral  particle  transport  methods  are  often  used  to  solve  complex  problems  in 
nuclear  science  and  engineering.  Several  Air  Force  programs  routinely  require  detail,  d 
numerical  models  of  neutron  and  photon  transport  in  nuclear  systems,  including  space 
nuclear  power  reactors  and  nuclear  effects  experiments.  The  innovation  of  new,  more 
efficient  numerical  schemes  for  solving  radiation  transport  problems  can  be  beneficial, 
often  providing  very  accurate  results  while  consuming  fewer  resources  than  required  by 
conventional  numerical  transport  methods.  Accurate  modeling  of  neutrons  and  gamma 
rays  in  a  given  system  can  be  achieved  using  a  discrete  ordinates  solution  to  the 
Boltzmann  transport  equation,  which  is  given  below  in  the  general  integro-differential 
form  for  neutrons  (Lewis  and  Miller,  1984:34-39): 

+  Cl  •  Vy(r,  Cl,E,t)  +  a(r,E)\\t(r,  Cl,E,  t)  = 

v  at 

s„,(r,  Cl,E,t)  +  JdE'J dCl'as(r, E’  ->  E,  Cl'  •  Cl,  )\| f(r,  Cl’,  E’,  t)  x  1 ' 

+X(E)jdE'vof(7,E'M7,E',t) 

where 

v  =  neutron  speed 

\\i  =  angular  neutron  flux 
Cl  =  unit  vector  in  direction  of  particle  motion 
r  =  coordinate  location  in  space 
E  =  neutron  energy 
t  =  time 

a  =  total  macroscopic  cross  section 
s„,  =  external  neutron  source 

ot  =  scattering  cross  section  causing  neutrons  to  arrive  in  E  and  Cl 


%  =  fission  neutron  energy  distribution  function 
vof  =  fission  neutrons  produced 
<j>  =  scalar  neutron  flux 


Functional  dependence  of  each  variable  is  defined  in  the  Boltzmann  transport  equation 
above.  The  right  hand  side  of  the  transport  equation  is  often  referred  to  as  the  source 
function  S(r,Cl,E,t),  which  incorporates  the  external,  scattering,  and  fission  source 
terms  implicitly.  A  discrete  ordinates  solution  to  equation  (1)  consists  of  evaluating  this 
equation  in  a  number  of  distinct  angular  directions  over  a  defined  spatial  mesh  over  a 
range  of  defined  energy  groups.  Suitable  angular  and  spatial  quadratures  are  used  to 
approximate  integral  particle  fluxes,  currents,  and  moments  (Lewis  and  Miller, 
1984:116-119). 

This  thesis  introduces  a  new  discrete  ordinates  spatial  quadrature  method  to  solve  the 
Boltzmann  transport  equation  in  slab  geometry.  This  approach,  the  "exponential 
characteristic"  (EC)  method,  assumes  an  exponential  form  for  the  position  dependence  of 
the  source  function,  defined  over  each  cell  in  a  way  that  preserves  the  zeroth  and  first 
spatial  moments  of  the  source  over  the  cell.  Cell  edge  fluxes  and  cell  flux  moments  are 
obtained  by  using  this  exponential  source  representation  in  a  characteristic  integration  of 
the  Boltzmann  transport  equation. 

In  addition  to  having  the  qualities  of  positivity,  continuity  over  a  cell,  and 
conservation  of  zeroth  and  first  spatial  source  moments,  numerical  testing  has 
demonstrated  that  this  new  exponential  characteristic  scheme  has  other  desirable 
properties: 
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•  Fourth  order  truncation  error  as  oAx  —>  0 

•  Proper  solutions  to  diffusive  problems,  as  tested 

•  Highly  accurate  over  very  coarse  meshes 


Ultimately,  an  efficient  method  that  proves  to  be  accurate  in  one  dimension  offers  the 
most  utility  if  easily  adaptable  to  multidimensional  coordinate  systems  so  that  more 
difficult  problems  might  be  solved.  Based  on  applications  of  the  exponential 
characteristic  method  in  one  dimension  with  isotropic  scattering,  presented  here,  this 
method  should  be  generalizable  to  multi-dimensional  geomr  tries. 
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In  Section  II,  necessary  background  information  is  discussed.  The  research  problem 
is  formally  stated  in  Section  III,  with  all  applicable  assumptions  and  limitations. 
Theoretical  development  of  the  exponential  characteristic  method  is  presented  in  Section 
IV,  including  methods  for  converting  equations  into  stable  computable  forms.  This  is 
followed  by  computer  program  development  and  code  validation  in  Section  V.  In 
Section  VI,  test  results  from  exponential  characteristic  solutions  to  three  test  cases  are 
evaluated  and  compared  to  other  discrete  ordinates  schemes.  Finally,  conclusions  and 
recommendations  are  given  in  Section  VII,  followed  by  several  appendices. 
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The  transport  equation  can  be  solved  only  numerically  for  complex  real  world 
applications,  e.g.,  reactor  core  and  shielding  problems;  therefore,  it  is  essential  that  any 
computer  algorithm  applied  be  efficient  and  accurate  to  a  degree  that  affords  a  solution 
within  a  reasonable  amount  of  time  and  cost.  There  are,  in  principle,  two  approaches  to 
spatial  quadrature  discrete  ordinates  schemes  used  in  solving  the  Boltzmann  transport 
equation: 

1)  Characteristic  schemes  -  based  on  the  known  streaming  properties  of  the 
transport  equation  (e.g.,  an  assumed  form  of  the  source  function) 

2)  Curve  fitting  schemes  -  based  on  an  assumed  form  of  the  flux  distribution 
function  \\i(r) 

Regardless  of  which  approach  is  used,  five  desirable  properties  of  spatial  quadrature 
schemes  to  considered  in  solving  the  Boltzmann  transport  equation  are: 

1)  accuracy  -  where  there  is  a  small  truncation  error 

2)  simplicity  -  where  a  small  number  of  numerical  operations  with  unknowns  from 
a  single  mesh  cell  are  involved 

3)  positivity  -  where  positive  flux  values  result  from  positive  sources  and  positive 
entering  boundary  fluxes 

4)  particle  conservation  -  where  particles  are  neither  created  nor  destroyed 
arbitrarily,  with  moment  balances  (e.g.  zeroth  moment,  first  moment,  etc.)  satisfied 
up  to  a  given  order 

5)  generalizabilitv  -  where  the  method  can  be  applied  to  other  geometries: 
Cartesian,  cylindrical,  or  spherical  (Lathrop,  1969:475-477) 


An  additional  quality  for  any  scheme,  introduced  by  Larsen,  requires  that  for 
problems  that  are  diffusive  (absorptions  and  gradients  of  the  angular  flux  are  small),  the 
solutions  resulting  from  any  given  transport  method  must  agree  with  solutions  arrived  at 
using  the  neutron  diffusion  equation.  This  criteria  is  defined  as  the  diffusion  limit.  For 
some  schemes,  satisfying  the  diffusion  limit  is  a  problem,  especially  when  the  spatial 
mesh  is  optically  thick  (Larsen,  1983:90-98).  It  is  difficult  for  a  numerical  scheme  to 
meet  all  of  the  above  criteria.  For  most  schemes,  maintaining  positive  fluxes  throughout 
a  problem  is  difficult  unless  optically  thin  cells  are  specified,  which  can  result  in  very 
high  computational  costs  (particularly  in  two  or  three  space  dimensions). 

A,  Diamond  Difference  Method 

Positivity  is  sacrificed  for  simplicity  in  the  implementation  of  the  popular  diamond 
difference  method.  Consider  a  one  dimensional  slab  in  the  x  direction  made  up  of  i  cells, 
where  any  ith  cell  has  length  Ac,  (between  x,  _ ,  and  x,)  with  homogeneous  material 
properties.  Diamond  difference  is  a  curve  fitting  scheme  where  one  assumes  \|/(.y),  is  a 
linear  function  across  a  cell  (Lathrop,  1969:482-484).  In  the  case  of  steady-state, 
mono-energetic,  non-multiplying,  isotropic  scattering  and  sources,  Cartesian  (slab) 
geometry,  the  Boltzmann  transport  equation  simplifies  to  (suppressing  the  subscript  i,  and 
noting  implicitly  that  each  angular  flux  ym(x)  is  along  a  specific  direction  cosine  (im): 

dyjx) 

P-  dx—  +  omy(x)  =  s(x)  (2) 

where  the  source  function  is 

sfO  =  ca<K.r) +  $„,(*)•  <3) 

The  scattering  ratio  in  equation  (3)  is  defined  as 
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(4) 


=  o, 

G 

The  scalar  flux  <j>(x)  for  any  spatial  quadrature  method,  defined  as  the  angular  flux 
integrated  over  all  angles,  is  given  in  slab  geometry  over  all  direction  cosines  ji  between 
-1  and  1  (corresponding  to  all  angles  between  it  and  0  radians,  respectively)  by 


4>(^)  =  Ji  V(*,4)y 


(5) 


Using  the  discrete  ordinates  method,  the  scalar  flux  is  approximated  by 


1  M 

<)>(*)  I  "mV/mW  (6) 

L  m  =  1 

where  M  angular  fluxes  are  evaluated  at  M  discrete  ordinates  p.m.  Each  direction  cosine 

is  from  a  quadrature  set  The  eight  direction  cosines  from  an  eight-point  single 

range  (S8)  Gauss-Legendre  quadrature  set,  used  for  all  computations  in  this  thesis,  are 
provided  in  Appendix  A. 

Equation  (2)  can  be  integrated  over  a  single  spatial  (mesh)  cell,  with  this  cell  defined 
locally  between  0  and  Ax.  When  jj.  >  0,  then  particles  flow  from  left  to  right,  and 
vji(0)  =  \\fL  at  x  =  0  and  \|/(Ar)  =  \|/s  at  x  =  Ar  (Note  that  L=Left,  R=Right,  and 
A=Average).  This  results  in  the  zeroth  moment  balance  equation: 

"  Vl)  +  =  sa A*  (?) 

where  the  cell  average  source  function  and  flux  are 


and 


(8) 
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(9) 


A  single  spatial  cell  is  shown  in  Figure  1. 


Figure  1.  A  Single  Spatial  Cell 


The  diamond  difference  method  approximates  equation  (9),  the  cell  average  angular 
flux  \\fA  (called  the  zeroth  moment  of  the  angular  flux)  by  the  arithmetic  average  of  the 
left  and  right  boundary  fluxes, 

(10) 

This  is  the  fundamental  assumption  of  the  diamond  difference  method.  Since  \\fL  and  SA 

are  either  given  from  an  initial  guess  or  are  available  from  a  previous  iteration,  the  basic 
diamond  difference  equation  is  obtained  by  combining  equations  (7)  and  (10): 


where  the  angular  optical  thickness,  in  mean  free  paths,  is  defined  as 


(11) 
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Ideally,  the  diamond  difference  method  yields  reasonable  results  and  second  order 
convergence  in  the  limit  of  optically  thin  cells,  where  oAx  — »  0  (Lewis  and  Miller, 
1984:132).  However,  even  if  \yL  and  SA  are  both  positive,  a  negative  angular  flux  is 
possible  on  the  right  cell  boundary  if  the  optical  thickness  is  greater  than  two.  Negative 
fluxes  can  occur  even  when  the  cell  thickness  Ax  is  small,  because  the  cross  section  a  can 
be  large,  and  the  direction  cosine  |i  can  be  small.  In  reality,  negative  fluxes  have  no 
physical  meaning.  In  addition  to  the  doubt  negative  fluxes  cast  on  the  usefulness  of 
results,  they  can  propagate  unfavorably  through  a  problem. 

This  can  be  treated  by  using  a  negative  flux  fixup,  where  if  <  0,  then  it  is  set  to 

zero.  The  zeroth  moment  balance,  equation  (7),  is  then  solved  for  \\fA.  While  practical, 
these  fixups  can  cause  undamped  oscillations  among  successive  iterations,  interfering 
with  or  even  preventing  convergence  (Alcouffe,  et  al.,  1979:1 1 1-1 15).  In  addition, 
negative  flux  fixups  increase  the  truncation  error  of  results,  can  cause  failures  in  diffusion 
acceleration  procedures,  and  can  ultimately  lead  to  incorrect  solutions  (Lathrop, 
1969:475). 

B.  Other  Spatial  Quadrature  Schemes 

Following  a  wide  use  of  the  diamond  difference  method  to  solve  transport  problems 
(in  spite  of  its  inherent  limitations),  many  researchers  sought  to  develop  alternative 
spatial  quadrature  schemes.  In  general,  other  methods  that  try  to  overcome  the 
difficulties  of  diamond  difference  are  less  straight  forward.  Nevertheless,  while  some  of 
these  schemes  are  more  intricate  to  encode,  many  achieve  more  accurate  results  than 


diamond  difference  over  a  coarser  spatial  mesh.  Therefore,  although  some  alternative 
methods  require  increased  computational  effort  per  cell  compared  to  diamond  difference, 
the  use  of  a  coarser  mesh  usually  lowers  overall  computation  times  and  memory  storage 
requirements  (This  is  particularly  important  in  two  and  three  dimensional  geometries). 
Some  alternatives  to  diamond  difference  are  listed  in  Table  1. 


Table  1 

Alternative  Spatial  Quadrature  Schemes  in  Discrete  Ordinates 


1.  (SC)  Step  Characteristic  (Lathrop,  1969) 

2.  (LD)  Linear  Discontinuous  (Hill,  1975) 

3.  (EM)  Exponential  Method  (Barbucci,  et  al.,  1977) 

4.  (LC)  Linear  Characteristic*  (Alcouffe,  et  al.,  1979) 

5.  (LN)  Linear  Nodal+  (Walters,  et  al.,  1981) 

6.  (LNA)  Augmented  Linear  Nodal  (Walters,  1986) 

7.  (SGF)  Spectral  Green’s  Function  (DeBarros,  et  al.,  1990) 

8.  (SA)  Step  Adaptive  (Mathews,  1990) 

9.  (LA)  Linear  Adaptive  (Mathews,  1990) 

10.  (MB)  Multiple  Balance  (Morel,  et  al.,  1990) 

+Linear  Nodal  is  identical  to  Linear  Characteristic  in  slab  geometry 


Although  most  are  improvements  to  diamond  difference,  each  method  in  Table  1  has 
distinct  advantages  and  disadvantages.  Detailed  descriptions  of  each  scheme  can  be 
found  in  the  literature.  Some  of  the  schemes  listed  in  Table  1  are  briefly  mentioned  here; 
others  are  left  for  discussion  in  later  sections  and  will  aid  in  the  development  of  the 
exponential  characteristic  method. 


The  linear  discontinuous  (LD)  method  was  developed  in  an  attempt  to  diminish  the 
number  of  negative  flux  fixups  required  by  the  diamond  difference  method.  The 
fundamental  assumption  used  in  the  LD  method  is  that  y  is  piecewise  linear  in  x,  but  is 
discontinuous  on  each  left  (entering)  cell  boundary  for  p  >  0,  and  discontinuous  on  each 
right  (entering)  cell  boundary  if  |i  <  0.  Similarly,  the  source  is  also  taken  to  be  piecewise 
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linear.  These  piecewise  representations  for  the  angular  flux  and  source  terms  are 
introduced  into  zeroth  and  first  spatial  moments  of  the  transport  equation.  This  generates 
two  equations  solvable  for  two  unknowns,  the  unknowns  being  the  two  discontinuous 
values  of  angular  fluxes  at  a  cell  edge.  The  LD  method  is  considered  to  be  a  near 
positive  method,  producing  fewer  negative  fluxes  than  the  diamond  difference  scheme  for 
a  given  problem.  Although  negative  fluxes  can  result,  the  linear  discontinuous  scheme 
does  not  invoke  a  fixup  and  achieves  second  order  convergence  in  the  limit  of  thin  cells. 
The  linear  discontinuous  scheme  is  algebraically  equivalent  to  the  linear  nodal  scheme  if 
a  Pade’  (1,2)  approximation  is  used  for  the  exponential  function  in  linear  nodal.  Because 
the  LD  scheme  continues  to  be  widely  implemented  in  transport  codes,  it  is  used  to  aid  in 
evaluating  the  performance  of  the  exponential  characteristic  method  (Hill,  1975), 
(Alcouffe,  et  al„  1979:1 17-1 18),  and  (Walters,  1981:1 15). 

Of  particular  note  is  the  exponential  method  (EM),  developed  by  Barbucci  and 
DiPasquantonio  in  1977  (not  to  be  confused  with  the  exponential  characteristic  scheme 
developed  in  this  effort).  Recall  that  the  diamond  difference  method,  referring  to 
equation  (10),  assumes  the  cell  average  angular  flux  is  the  arithmetic  mean  of  the 
entering  and  exiting  fluxes.  Similarly,  the  exponential  method,  also  a  curve-fitting 
scheme,  assumes  the  cell  average  flux  is  the  geometric  mean  of  entering  and  exiting 
fluxes,  which  results  in  an  exponential  functional  form  for  all  angular  fluxes.  This 
method,  while  positive,  has  fundamental  difficulties  in  source  regions  and  near  vacuum 
boundaries.  In  addition,  exact  solutions  are  always  overestimated,  even  in  the  limit  of 
thin  cells  (Barbucci,  et  al.,  1977).  Because  of  these  difficulties,  no  comparisons  with  the 
exponential  method  are  made. 


In  slab  geometry,  there  are  no  real  differences  between  the  linear  characteristic 
method  (LC)  (discussed  in  a  later  section)  and  the  linear  nodal  (LN)  method,  except  in 
the  implementation  of  a  fixup  (the  reader  is  referred  to  the  literature  for  treatment  in  other 
geometries).  The  LN  method  actuates  a  fixup  if  the  magnitude  of  the  first  moment  of  the 
scalar  flux  exceeds  that  of  the  scalar  flux  zeroth  moment,  while  the  LC  method 
implements  a  fixup  based  on  the  magnitudes  of  the  source  moments.  In  the  limit  of  thin 
cells,  LC  and  LN  are  equivalent  in  slab  geometry.  The  augmented  linear  nodal  scheme 
(LNA)  implements  the  diamond  difference  relation  for  specific  terms  in  two  and 
three-dimensional  moment  balance  equations  to  increase  computational  efficiency  (with  a 
small  penalty  in  accuracy).  In  any  case,  the  LNA  method  is  not  always  positive,  and  a 
fixup  is  still  required  (Walters,  et  al,  1981),  (Alcouffe,  et  al.,  1979),  and  (Walters,  1986). 
Because  LN  is  equivalent  to  LC  in  slab  geometry,  no  other  reference  to  it  is  made. 
Furthermore,  as  the  LNA  scheme  is  neither  absolutely  positive  nor  as  accurate  as  the 
original  LN  scheme,  it  is  not  discussed  further. 

The  spectral  Green’s  function  method  (SGF)  has  no  spatial  truncation  error  in  slab 
geometry.  In  addition  to  using  the  standard  transport  balance  equations,  it  implements  a 
non-standard  Green’s  function  for  \j/A,  which  affords  an  exact  solution  in  slab  geometry. 

In  slab  geometry,  this  method  requires  significantly  more  computer  memory  than  most 
other  schemes.  Furthermore,  the  method  does  not  appear  to  be  extensible  to  other 
geometries  (De  Barros,  et  al.,  1990).  For  these  reasons,  the  SGF  method  is  not 
considered  for  comparison  to  the  EC  scheme. 

Multiple  balance  (MB)  is  a  curve  fitting  scheme  that  uses  the  Boltzmann  transport 
equation  as  the  principle  balance  equation  and  auxiliary  balance  equations  restricted  to 
coupled  half-cell  domains.  The  MB  method  employs  a  Pade’(0,2)  approximation  to  the 
exponential  function.  While  this  scheme  has  many  desirable  properties,  it  requires 
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computations  over  cell  sub-domains,  is  limited  to  second  order,  and  is  not  always  positive 
when  scattering  is  present  in  a  given  region  (Morel,  et  al.,  1990).  Therefore,  no 
performance  comparison  between  multiple  balance  and  exponential  characteristic  is 
made. 

C.  Characteristic  Schemes 

i.  Step  Characteristic  Method 

The  step  characteristic  (SC)  method  uses  a  characteristic  integration  of  the  transport 
equation,  which  assumes  the  source  function  takes  the  form  of  a  step  value  s(x)  =  SA  in 
each  cell.  Thus,  analytically  integrating  equation  (2)  along  each  characteristic  direction 
(constant  |i)  yields  solutions  for  the  flux  distribution  tuid  zeroth  moment.  This  method  is 
simple,  positive,  and  conserves  zeroth  moments,  but  is  difficult  to  generalize  in 
multidimensional  geometries.  Also,  this  method  demonstrates  less  than  second  order 
truncation  error  as  cells  become  optically  thin  (oAx  -»  0)  (Lathrop,  1969:475). 

ii.  Linear  Characteristic  Method 

The  linear  characteristic  method  is  a  logical  extension  of  the  step  characteristic 
method.  Instead  of  approximating  the  source  function  as  a  step  value  throughout  a  cell, 
this  method  uses  a  linear  approximation  of  the  source  function  in  a  characteristic 
integration  of  the  transport  equation: 

s(x)  =  SAP0(x)  +  SrP](x)  (13) 
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where  the  zeroth  and  first  moments  of  the  source  function  are  SA  and  Sx,  respectively, 

and  P0  and  Pt  are  Legendre  polynomials.  These  moments  are  initially  obtained  by 
integrating  the  product  of  true  source  function  and  the  zeroth  and  first  order  Legendre 
polynomials.  The  true  source  function,  as  given  in  equation  (3),  is: 

s(*)  =  ca<KO +  *„.(*)  (3) 

The  Legendre  polynomials  P0  and  P„  defined  so  as  to  be  orthogonal  over  the  interval 
0  <x  <  Ax  ,  are: 


W  =  1 

04) 

(  2x 

Pl(x)=  — -1 

Integration  is  carried  out  using  equations  (3),  (14),  and  (15)  according  to 

(15) 

Sa=ZxJo  s(x)F°(x)dx 

(16) 

(17) 

These  result  in 

SA  =  co<pA  +SalA 

(18) 

Sx  =  c  otyx  + 

(19) 
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Now  substituting  the  linear  approximation  for  s (x)  from  equation  (13)  into  the 
transport  equation  (equation  (2))  and  integrating  analytically  yields  a  solution  for  the 
angular  flux  distribution  in  the  cell  (where  y(0)  =  \j/L  is  the  left  entering  flux  for  a 
specified  p  >  0): 


(20) 


Evaluating  this  equation  at  x  =  Ax  yields  the  right  boundary  angular  flux  \\iR.  Integration 

of  equation  (20)  can  be  performed  to  obtain  the  zeroth  and  first  angular  flux  moments 
using 


1  f** 

Va=^Jo  V(x)P0(x)dx 
3 


(21) 


(22) 


The  expressions  for  and  4>x  in  the  zeroth  and  first  source  moments  (equations  (18)  and 

(19)),  are  actually  the  angular  flux  moments  integrated  over  all  angles  (in  slab  geometry) 
according  to: 


♦a  =  Ji  V,(M)y 

(23) 

(24) 
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Naturally  all  angular  flux  integrations  are  performed  numerically  using  the  selected 
angular  quadrature  set  under  the  discrete  ordinates  approximation  (see  equations  (5)  and 
(6)). 


Figure  2.  Iteration  Sequence  in  Linear  Characteristic  Spatial  Quadrature 


A  schematic  of  the  progression  from  iteration  to  iteration  in  linear  characteristic 
spatial  quadrature  is  given  in  Figure  2.  Since  the  source  moments  are  either  known  from 
an  initial  guess  or  can  be  evaluated  from  a  previous  iteration,  these  yield  new  fluxes  and 
flux  moments,  which  provide  new  source  moments  for  each  cell.  This  process  continues 
until  convergence  is  achieved  to  a  specified  tolerance.  Logically,  this  process  is  often 
referred  to  as  fixed  point  iteration  of  the  source- scattering  term. 
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If  |  SJ  >  SA :  which  can  occur  when  cells  are  optically  thick,  the  linear  source 

distribution  equation  (13)  is  negative  at  one  end.  To  preserve  positivity,  a  negative 
source  fixup,  also  known  as  a  rotational  fixup,  is  required.  This  source  fixup  restricts  the 
absolute  value  of  Sx  so  that  it  is  never  larger  than  the  zeroth  source  moment  SA.  While 
this  fixup  insures  positive  sources  (and  hence  positive  fluxes),  it  violates  first  spatial 
moment  conservation  of  the  Boltzmann  transport  equation  (and  forces  truncation  error  to 
increase).  This,  in  turn,  can  cause  numerical  diffusion;  particles  are  spatially  shifted  and 
can  appear  where  they  should  not.  The  extent  to  which  numerical  diff  usion  occurs  in  LC 
depends  on  the  impact  of  limiting  the  magnitude  of  Sx.  In  short,  this  source  fixup  can 
introduce  problems  analogous  to  those  caused  by  negative  flux  fixups  in  the  diamond 
difference  method. 

If  the  rotational  fixup  is  never  used,  the  linear  characteristic  method  globally 
conserves  the  zeroth  and  first  spatial  moments  of  the  Boltzmann  transport  equation.  The 
zeroth  spatial  moment  balance  equation  in  slab  geometry,  the  same  balance  equation  used 
in  the  diamond  difference  approximation  (equation  (7)),  is 

H(V*  -  Vi)  +  oAx-Va  =  SA  At  (25) 

The  first  spatial  moment  balance  equation  is 

3|i(\j/t  -  2  \rA  +yR)  +  cAr  y,  =  Sx  At  (26) 

These  moment  balances  of  the  Boltzmann  transport  equation  are  easily  derived  by 
integrating  equation  (2)  multiplied  by  the  properly  normalized  zeroth  and  first  order 
Legendre  functions,  equations  (14)  and  (15),  between  the  limits  0  and  At.  Over  coarse 
meshes,  the  linear  characteristic  method  is  second  order.  In  the  limit  of  vanishingly  thin 
cells,  truncation  of  the  LC  method  is  at  least  third  order  (Alcouffe,  et  al.,  1979:1 1 1-127). 
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In  spite  of  potentially  requiring  a  fixup  in  some  cells  (depending  upon  the  problem), 
the  linear  characteristic  method  has  produced  better  results  than  the  diamond  difference, 
linear  discontinuous,  exponential,  and  step  characteristic  methods  in  several  test  cases 
(Alcouffe,  et  al.,  1979:126-127).  In  addition,  it  satisfies  the  correct  diffusion  limit 
(Larsen,  1983:90).  Therefore,  LC  remains  a  worthy  scheme  with  which  to  compare. 

D.  Adaptive  Schemes 

Mathews  recently  developed  two  methods,  step  adaptive  (SA)  and  linear  adaptive 
(LA)  (Mathews,  1990:419-457).  These  are  logical  improvements  to  the  step  and  linear 
characteristic  methods,  respectively.  These  new  adaptive  methods  differ  from  earlier 
schemes  in  that  source  functions  are  strictly  positive  throughout  the  cell  and  require  no 
rotational  fixup.  This  is  accomplished  by  partitioning  the  domain  of  each  cell  into 
smaller  sub-domains.  The  source  function  is  represented  by  different  functions  over  each 
sub-domain.  The  operative  "adaptive"  for  these  schemes  refers  to  the  precise  partitioning 
of  the  sub-domains  in  a  way  that  forces  the  source  representation  to  conserve  0th  and  1st 
moments  in  a  spatial  cell.  Hence,  the  source  representation  is  determined  by  a  moments 
matching  process.  Due  to  the  intrinsic  positivity  of  this  method,  no  fixups  are  required, 
and  zeroth  and  first  moments  of  sources  and  flux  distributions  are  always  conserved,  as 
are  the  balance  equations  (equations  (25)  and  (26)).  Note  that  these  methods,  while 
positive,  are  non-linear.  This  non-linearity  has  not  been  observed  to  interfere  with 
convergence  of  the  scattering  iteration. 

Not  surprisingly,  Mathews  found  SA  and  LA  to  be  in  excellent  agreement  with 
conventional  methods  in  test  problems  with  thin  cells.  However,  when  thick  cells  were 
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used,  these  adaptive  schemes  formidably  outperformed  the  other  methods.  Truncation 
errors  of  third  and  fourth  order  are  realized  in  the  limit  of  very  thin  cells  with  SA  and  LA, 
respectively. 

i.  Step  Adaptive  Method 

In  the  step  adaptive  method,  a  step  function  approximates  the  source  function  over  a 
portion  of  the  domain,  and  is  truncated  to  zero  in  the  remainder.  The  location  of  the 
truncation  is  dependent  on  the  magnitudes  of  the  zeroth  and  first  moments  and  by  the 
sign  of  the  first  moment. 

For  example,  if  Sx  >  0,  then 


where 


sC 0  =  0 

5, 


sfr)  = 


(1-P) 


0<x<pAr  (27) 
pAr  <  x  <  Ax  (28) 


If  Sx  <  0,  then 


0  <  p  <  1  (29) 


S* 

s(x)  =  ■  -  0  <*  <  (1  -p)Av  (30) 

$00  =  0  (1  -p)Ar  <x  <  Ax  (31) 
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ii.  Linear  Adaptive  Method 


The  linear  adaptive  method  works  in  a  manner  similar  to  SA,  but  is  more  accurate.  If 
|  Sx |  <  SA,  then  s(.r)  takes  the  form  of  equation  (13),  as  in  the  LC  method.  However,  if 
Sx  >  SA ,  then  $(*)  is  set  to  zero  from  the  left  edge  of  the  cell,  up  to  a  distance  (1  -  x)Ar, 
where  x  =  3(1  -  p)/2.  From  this  point,  it  rises  linearly  to  a  maximum  at  the  right  edge.  A 
mirror  image  of  this  is  used  if  Sx  <  -(S^).  This  insures  the  source  function  is  continuous, 
positive,  and  piecewise  linear  in  the  cell;  again,  zeroth  and  first  moments  are  conserved. 
In  the  limit  of  vanishingly  thin  cells  in  slab  geometry,  the  linear  adaptive  method  is 
equivalent  to  the  linear  characteristic  method  (Mathews,  1990:419-457). 

E.  Merits  of  an  Exponential  Characteristic  Method 

The  true  power  of  the  step  and  linear  adaptive  schemes  is  in  their  ability  to  retain 
accuracy,  positivity,  and  conservation  with  coarse  meshes.  This  suggests  the  possibility 
of  solving  complex  transport  problems  once  thought  too  large  or  costly  for  conventional 
computer  systems.  One  potential  difficulty  in  using  SA  and  LA  involves  adapting  the 
source  function  to  match  moments  over  a  sub-domain.  If  posed  with  a  problem  with 
strong  and  weak  fluxes  appearing  at  opposite  boundaries  of  a  given  cell  and  |  Sx\  >  SA,  the 
adapted  source  function  is  preferentially  biased  to  de-emphasize  the  weak  boundary  flux, 
while  the  strong  boundary  flux  is  over-emphasized.  An  additional  obstacle  in  SA  and  LA 
is  simply  the  added  complexity  of  using  a  piecewise  source  representation  in 
sub-domains  across  a  mesh  cell.  Depending  on  the  problem  being  solved,  use  of 
sub-domains  can  require  more  floating  point  operations  than  might  be  necessary  for  a 
new  method  that  uses  a  continuous  source  representation  across  an  entire  cell.  (Use  of 
sub-domains  are  also  more  difficult  to  treat  in  two  or  three  space  dimensions). 
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Still,  any  newly  proposed  cell-continuous  source  function  must  provide  results  that 
retain  the  features  of  accuracy,  positivity,  and  conservation  of  at  least  zeroth  and  first 
moments  (even  over  a  coarse  mesh).  One  source  function  that  satisfies  all  these 
requirements  is  the  exponential  function 

s(jc)  =  aexp[b;c]  (32) 

where  a  and  b  are  constants  to  be  chosen  so  as  to  conserve  the  zeroth  and  first  spatial 
moments  (i.e.  to  match  the  specified  moments,  SA  and  Sx,  from  the  outer  iteration).  This 
scheme  is  called  the  exponential  characteristic  (EC)  method.  While  the  constants  in  the 
source  function  in  equation  (32)  are  defined  by  moments  matching  similar  to  that  used  for 
S  A  and  LA,  the  exponential  characteristic  method  is  not  an  adaptive  scheme;  it  will  not 
use  sub-domains  across  a  spatial  cell.  An  added  feature  of  this  source  function  is  its 
similarity  to  panicle  distributions  in  strongly  absorbing  media.  Like  the  SA  and  LA 
schemes,  the  EC  method  is  non-linear  (as  in  schemes  which  use  a  fixup),  and  therefore 
does  not  preserve  the  linearity  of  the  transport  equation  in  the  sense  of  the  super-position 
of  solutions. 
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III.  Problem 


The  problem  addressed  here  is  the  theoretical  development,  computer  program 
implementation,  validation,  testing,  and  evaluation  of  a  discrete  ordinates  exponential 
characteristic  spatial  quadrature  scheme  to  numerically  solve  the  Boltzmann  transport 
equation.  The  EC  method  is  evaluated  based  on  a  comparison  of  numerical  results  to 
results  from  other  spatial  quadrature  schemes. 


A.  Scope 

The  "exponential  characteristic"  method  is  derived  using  a  characteristic  integration 
of  the  Boltzmann  transport  equation  with  an  exponential  function  as  the  assumed  form  of 
the  source  distribution.  Although  the  method  derived  here  could  be  implemented  more 
generally,  for  simplicity  of  programming,  the  implementation  and  testing  performed  here 
is  restricted  to  the  conditions  stated  in  Table  2. 


Table  2 

Conditions  Used  in  Testing  the  EC  Method 


1.  Steady  State 

2.  Mono-Energetic 

3.  Non-Multiplying 

4.  Isotropic  Scattering  and  Sources 
in  the  Laboratory  System 

5.  Slab  Geometry 


Validation  and  testing  of  the  EC  scheme  is  performed  using  test  cases  with  multiple 
regions  containing  combinations  of  scattering  and  absorbing  media.  In  each  case, 
conservation  of  zeroth  and  first  spatial  moments  is  numerically  verified  using  the  balance 
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equations  (equations  (25)  and  (26)).  In  addition,  an  overall  scalar  balance  equation  (the 
steady  state  Boltzmann  equation  integrated  over  all  angles)  is  also  numerically  verified. 
This  overall  balance  equation  is,  in  slab  geometry, 

7,  (Ax) -7,(0)  +  oAx<t>A  =  SaAx  (33) 

where 

7,(Ax)  =  net  current  at  right  cell  boundary 

7,(0)  =  net  current  at  left  cell  boundary 
a  =  total  macroscopic  cross  section 
Ax  =  measured  cell  width 
<t>A  =  scalar  flux  zeroth  moment 
SA  =  source  function  zeroth  moment 


Note  that  the  boundary  currents  7, (Ax)  and  7,(0)  are  given  in  slab  geometry  by 


and 


,(Ax)=  f  |i\|/(Ax,n)~  I  limH-m\j/m(Ax) 

J-\  Z  m  =  1 


(34) 


7,(0)=  f  qv(0,q)^=  I  )tmH’nym(0)  (35) 

J- 1  Z  m  -  1 

As  before,  all  angular  integrations  are  carried  out  using  the  discrete  ordinates 
approximation. 
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Evaluation  of  the  performance  of  the  exponential  characteristic  method  is,  for  any 
given  test  problem,  achieved  by  comparing  numerical  results  obtained  from  the  EC 
method  with  results  obtained  from  conventional  spatial  quadrature  schemes  in  slab 
geometry.  Comparisons  are  made  among  the  schemes  listed  in  Table  3. 


Table  3 

Spatial  Quadrature  Schemes  Used  in  Test  Problems 


(EC)  Exponential  Characteristic 

(DD)  Diamond  Difference 

(DDF)  Diamond  Difference  with  Fixup 

(LD)  Linear  Discontinuous 

(LC)  Linear  Characteristic 

(SA)  Step  Adaptive 

(LA)  Linear  Adaptive 


I  A  relative  convergence  tolerance  of  10  5  using  Gauss-Legendre  angular  quadrature  with 

eight  direction  cosines  (S8)  is  used  to  solve  all  problems.  Only  a  numerical  comparison  of 
results  is  provided;  no  rigorous  mathematical  proof  of  the  efficiency  of  the  EC  scheme 
versus  other  schemes  is  made. 


Because  equations  (25),  (26),  and  (33)  are  numerically  verified  for  each  quadrature 
scheme,  including  exponential  characteristic,  a  check  on  all  angular  and  scalar  quantities 
is  available.  If  numerical  rounding  errors  or  catastrophic  cancellations  have  degraded 
solutions,  this  is  immediately  obvious  from  the  inspection  of  the  relative  maximum 
differences  in  these  balance  equations.  (Note  that  the  first  spatial  moment  balance 
equation  (26)  is  not  used  for  DD  or  DDF). 
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B.  Error  Norm  and  Error  Treatment 


For  easy  interpretation  and  evaluation  of  results,  the  error  norm  given  by  Mathews  as 


Err  (approx)  = 


|  approx  -  exact  \ 
(\approx\  - 1  exact]  )/2 


(36) 


is  used  throughout  all  studies  to  compare  all  cell  edge  fluxes,  currents,  and  region 
moments,  as  appropriate.  In  addition,  this  same  norm  is  used  to  compare  balance 
equations  (25),  (26),  and  (33)  for  each  quadrature  scheme  during  computer  execution, 
where  the  right  and  left  hand  sides  of  these  equations  are  taken  to  be  approx  and  exact , 
respectively.  While  many  authors  have  used  slightly  different  error  norms,  this  norm  is 
stable  for  zero  or  negative  results,  and  approaches  a  conventional  relative  error  estimate 
when  the  approximate  value  is  close  to  the  exact  value  (Mathews,  1990:446).  By 
inspection,  it  is  clear  that  the  maximum  value  produced  by  this  error  norm  is  two. 


Average  and  maximum  observed  errors  are  determined  using  error  values  computed 
from  equation  (36).  An  average  pointwise  relative  error  Q  provides  a  measure  of  the 
accuracy  of  all  computed  values  over  an  entire  problem,  and  is  given  by 


1  N 

Q=jj  jL  Err  (approx  )„  (37) 

where  N  values  for  all  fluxes,  currents,  etc,  are  computed  for  a  given  spatial  quadrature 

method.  The  maximum  observed  error  is  the  maximum  relative  error  observed  in  any 
quantity  for  a  given  solution  to  a  problem;  therefore,  it  provides  a  "worst  case" 
calculation  error  for  a  spatial  quadrature  method  used  to  solve  the  problem. 

As  for  the  "exact"  numerical  solution  from  which  errors  are  computed  in  equation 
(36),  this  is  the  reference  solution  obtained  using  the  linear  characteristic  method  with  a 
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very  fine  spatial  mesh  (where  further  refinement  of  the  mesh  yielded  no  discemable 
change  in  the  numerical  answer).  Unless  otherwise  stated,  a  relative  tolerance  of  10"5  is 
always  used  for  convergence. 


C.  Computational  Efficiency  and  Cost 

In  defining  computational  efficiency,  an  operational  definition  analogous  to  that 
stated  by  Alcouffe,  et  al.,  is  used.  This  holds  that  one  method  is  more  computationally 
efficient  than  another  if,  for  the  same  relative  cost,  the  first  method  produces  a  more 
accurate  solution  (Alcouffe,  et  al.,  1979:1 13).  Here,  the  relative  cost  in  solving  any  test 
case  is  determined  by  the  size  of  the  spatial  mesh  and  the  execution  time  required  for 
each  spatial  quadrature  method  to  reach  a  converged  solution.  Criteria  for  convergence  is 
a  comparison  of  scalar  fluxes  in  all  cells  computed  before  and  after  the  last  iteration.  If 
the  relative  change  (using  equation  (36))  in  these  scalar  fluxes  is  no  larger  than  1(T5,  then 
the  convergence  test  is  satisfied. 


The  average  error  equation  (37)  is  used  to  determine  a  mesh  size  ratio,  or  MSR.  The 
MSR  is  defined  as 


MSR= 


Cells  (Q)^' 
Cells  (Q)ec 


(38) 


where 

Cells  (Q  )xx  =  number  of  cells  required  to  achieve  an  average  error  Q  using 
spatial  quadrature  scheme  XX 

Celts  (Q)ec  ~  number  of  cells  required  to  achieve  an  average  error  Q  using  the 
EC  scheme 

XX  =  DD,  DDF,  LD,  LC,  SA,  or  LA  as  in  Table  3. 
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Another  useful  ratio  is  the  execution  time  ratio,  or  ETR.  This  is  defined  as 


'  Exec  Time  ( Cells(Q  ))yv ' 

ETR=  — — —  7  ^  (39) 

Exec  Time  (Cells(Q  ))EC 

where 

ExecTimexx  =  execution  time  required  for  method  XX  to  achieve  convergence 
using  Cells  (Q)xx 

ExecTimeEC  =  execution  time  required  for  the  EC  method  to  achieve 
convergence  using  Cells(Q)£C 

XX  =  DD,  DDF,  LD,  LC,  SA,  or  LA  as  in  Table  3. 

Memory  storage  requirements  are  proportional  to  the  number  of  mesh  cells  needed  to 
achieve  a  specified  error,  and  computational  effort  is  related  to  the  execution  time 
necessary  for  convergence  using  the  number  of  cells  designated.  Thus,  while  the  average 
pointwise  error  Q  provides  a  measure  of  overall  solution  accuracy,  the  MSR  and  ETR 
provide  a  comparison  of  memory  storage  requirements  and  computational  effort, 
respectively,  for  any  method  XX  against  the  memory  storage  and  computational  effort 
required  for  EC.  The  higher  the  value  of  either  the  MSR  or  ETR  above  unity,  the  greater 
the  benefit  from  using  the  exponential  characteristic  scheme.  A  mesh  cell  ratio  or 
execution  time  ratio  less  than  unity  for  a  given  method  indicates  fewer  cells  or  a  shorter 
execution  time,  respectively,  than  for  EC. 

D.  Source  Iterations 

Two  approaches  to  the  numerical  iteration  of  the  source  function  are  used  in  the 
study  of  each  test  case.  The  first  approach  is  a  fixed  point  iteration  of  the  source 
function,  where  scalar  flux  moments  are  updated  between  each  iteration  to  provide  zeroth 
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and  first  source  moments  (see  Section  Il.C.ii).  The  second  approach  is  the  method  of 
successive  scatters.  Iteration  using  successive  scatters  differs  from  fixed  point  iteration 
in  that  any  external  sources  and  incident  currents  present  are  used  to  inject  a  flood  of 
neutrons  into  the  system  during  the  first  iteration  only;  this  yields  the  first  flight  flux. 
Subsequent  iterations  deal  only  with  the  scattered  components  of  this  first  flight  flux. 
Initially,  the  source  term  is 


s(x)  =  sexl(x)  (40) 

For  the  second  iteration  and  thereafter,  the  source  term  is 

s(x)  =  cG§FIlghl(x)  (41) 

where  §Fllghl(x)  is  the  scalar  flux  generated  from  scattering  by  the  previous  flight.  As  the 

scattering  process  continues,  running  sums  of  all  fluxes  and  currents  are  accumulated 
until  the  scattered  particle  fluxes  become  negligible  and  convergence  is  achieved. 

The  treatment  of  the  source  function  in  successive  scatters  (given  by  equations  (40) 
and  (41))  differs  from  the  fixed  point  iteration  source  term,  which  always  includes 
components  from  both  external  sources  and  scatters,  as  in  equation  (3).  Since  all 
problems  treated  here  are  time  independent,  either  method  provides  essentially  the  same 
solution.  However,  the  method  of  successive  scatters  can  be  more  accurate  for  a  given 
mesh  size.  Because  the  external  source  term  often  dominates  the  source  function  in  fixed 
point  iteration,  contributions  from  scattered  radiation  can  be  masked;  this  does  not  occur 
in  successive  scatters,  as  only  scattering  from  a  previous  flight  contributes  to  the  source 
term  (and  thus  to  source  moments  for  each  iteration).  This  affects  the  accuracy  of  a 
solution,  to  an  extent  that  depends  on  the  spatial  quadrature  scheme  used.  This  is 
demonstrated  in  a  later  section. 
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The  exponential  characteristic  method  is  augmented  into  a  slab  geometry  discrete 
ordinates  computer  code  that  has  the  capability  of  executing  all  other  required  spatial 
quadrature  schemes.  Programs  were  written  in  Microsoft  QuickBASIC  4.5  for  execution 
on  an  IBM  or  compatible  personal  computer. 
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IV.  Theoretical  Development 

The  exponential  characteristic  method  is  derived  in  a  manner  very  similar  to  that  used 
for  the  linear  characteristic  method,  except  that  the  assumed  form  of  the  source  function 
is  an  exponential.  Again,  a  characteristic  integration  is  performed  over  a  single  spatial 
cell,  as  in  Figure  1. 

A.  Flux  Distributions 

Recalling  the  steady-state,  one  energy,  transport  equation  in  slab  geometry  (equation 

(2)): 


dym{x) 

dx—  +  omV(x)  =  s(x)  (2) 

Assume  the  isotropic  source  function  for  the  exponential  characteristic  method  has  the 
form  of  equation  (32): 

s(x) =  o  exp[Z?.v  ]  (32) 

Recall  that  a  and  b  are  to  be  determined  by  zeroth  and  first  source  moment  matching. 

Integrating  equation  (2)  by  integrating  factor  along  a  direction  cosine  |i,  and  noting  that  at 
the  left  hand  boundary,  \j/(0)  =  \| fL,  results  in 


V(jr)  =  Vtexrf 


'  —ox ' 

+  fVtexpT-^--0] 

.  4  . 

l  H  J 

dx_ 

M 


(42) 


Substituting  equation  (32)  into  equation  (42)  and  carrying  out  the  integration  yields  an 
expression  for  the  angular  flux  distribution  in  a  cell. 
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¥W  =  Vlexp[^]+(  ^  )(expM-exp[^])  (43) 

Since  \\iR  is  defined  as  y(Ar ),  evaluating  equation  (43)  at  x  =  Ax:  yields  the  angular  flux  at 
the  right  cell  boundary: 


\]/R=\\fLe\ 


'  -cAx ' 

l  “I 


b\i  +  G 


exp[6Ar] 


'  -cAx ' 
-exp  - 

L 


B.  Flux  Moments 

As  mentioned,  a  requirement  imposed  on  the  exponential  characteristic  scheme  is 
that  zeroth  and  first  moments  are  to  be  globally  conserved.  The  zeroth  angular  flux 
moment,  recalling  equation  (21),  is 

1  f* 

Va=^:Jo  V{x)P«{x)dx  (21) 

Substituting  equation  (43)  into  this  equation  and  integrating  yields  the  simplified  zeroth 
angular  flux  moment: 

T  -oAr'H  .  f  )/t 

M iA=—r~  1-exp  — —  +  —  (1-expfrAv]) 

1  \i  \)  (45) 

f  ~a\i  Y  r i  a  "l  r  -<yAx  Y 

{(b\x  +  o)oA.x 1  H  J  J 

Similarly,  making  use  of  equation  (22)  and  again  substituting  equation  (43)  yields  the 
angular  flux  first  moment  equation: 


31 


¥r  = 


6  a 


oAv  J 
3  a 


1  +  exd  ■ 


V 


'  -oAx ' 

\ 

[  6ViH2] 

r 

'  -a At  "j N 

>  - 

L  H  J 

+ 

S 

l  J 

1  -  exp 

V 

.  V-  J  J 

bo  Ax 


f 


(1  +  exp[6Av])  + 


—6(3(1 


V  b2oAx2  j 


\  ( 

(1  -exp[6Av])  + 


\b  c^Ar 
-3  a\i 


(1  -exp[6Ar]) 


X 


~6a\r 


(6|i  +  o)crAr 


^(dp  +  a)crAv 

V 

exp[6Av]  -expl 


exp[&Av]  -exd 


-aAr 


L  H  J) 


j 


—oAx 


(46) 


C.  Source  Moments 

Maintaining  the  requirement  that  source  moments  be  conserved,  the  zeroth  source 
moment  is  defined  as  given  in  equation  (16): 

S/4=iJo  s{x)P°{x)dx  (16) 

Substituting  equation  (32)  into  this  equation  yields  a  relationship  between  a,  b,  and  SA: 

5*=^(eXpl^Vl_1)  (47) 

Again,  the  first  source  moment  is  defined  as  in  equation  (17), 

s(x^P^x)dx  (17) 

and  substituting  equation  (32)  here  yields  a  relationship  between  a,  b<  and  Sx : 


5  = 


•,  l[Z;Ar(exp[£>Av]  +  1)  +  2(1  -exp[/?Ai])] 

Vb  Ar*; 


(48) 
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To  verify  the  expressions  for  the  moments  of  all  fluxes  and  sources,  equations  (44)  to 
(48)  were  substituted  into  the  balance  equations  (25)  and  (26)  to  yield  mathematical 
identities. 

D.  Method 

SA  and  Sx  are  both  known  either  as  an  initial  guess  or  as  evaluated  from  the  previous 

iteration.  The  equations  for  SA  are  as  derived  for  the  linear  characteristic  method. 
Recalling  equation  (18)  (where  c  is  given  in  equation  (4)): 

SA  =  co<bA  +  SalA  (18) 

the  zeroth  scalar  flux  moment,  using  equations  (21),  (23),  and  the  discrete  ordinates 
approximation,  is 


*A  I  (49) 

L  m  -  1 

The  external  zetoth  source  moment  is  also  defined  as  in  the  LC  method: 

1  C * 

s„,a  =^Jo  )/>(, )dx  (50) 

Note  that  the  value  specified  for  an  external  source  in  a  region  for  any  test  case  is 
assumed  to  be  the  zeroth  external  source  moment  (average  value  of  the  source  in  the 
region),  regardless  of  the  spatial  quadrature  method.  The  first  external  source  moment  is, 
for  the  problems  tested  here,  assumed  to  be  zero.  The  corresponding  equations  defining 
Sx  are,  recalling  equations  (19),  (22),  and  (24),  again  using  the  discrete  ordinates 
approximation 
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Sx  =  co<j>x  +  S„„ 

(19) 

1  « 

^  2  m  =  1  M’'n¥x'" 

(51) 

3  r * 

^Jo 

(52) 

The  Legendre  functions  are  as  in  equations  (14)  and  (15).  A  schematic  of  the  progression 
from  iteration  to  iteration  in  exponential  characteristic  spatial  quadrature  is  given  in 
Figure  3. 


34 


In  each  iteration,  SA  and  Sx  are  computed  from  equations  (18)  and  (19),  and 

coefficients  a  and  b  are  determined  by  root  solving  the  source  moment  relationships 
given  in  equations  (47)  and  (48).  Values  for  flux  moments  and  fluxes  are  then  computed 
using  equations  (44),  (45),  and  (46).  This  process  continues  over  all  directions  and  cells 
until  a  converged  solution  is  reached.  Since  the  source  function  is  never  negative,  no 
fixups  are  required,  and  the  EC  method  globally  conserves  zeroth  and  first  spatial 
moments.  Because  of  this,  the  balance  equations  given  in  equations  (25)  and  (26)  are 
satisfied  in  each  cell  for  each  discrete  ordinate.  In  view  of  this,  the  overall  balance, 
equation  (33),  must  also  be  satisfied.  This  is  true  unless  numerical  rounding  errors, 
catastrophic  cancellations  (overflows/underflows),  or  use  of  numerical  formulations  that 
are  ill-conditioned  compromise  the  solutions. 


E.  Conversion  to  Computable  Forms 

The  equations  derived  for  the  exponential  characteristic  method  are  not,  in  all  cases, 
numerically  stable.  For  example,  consider  equation  (44),  the  angular  flux  on  the  right 
boundary. 


vi/K=ytex 


P 


-oAlV 

~  j 


+ 


a 


b^i  +  c 


|^exp[6Ax]  -  ex 


(44) 


When  the  constant  b  is  zero  or  negative,  the  second  term  in  equation  (44)  can  become 

unstable.  By  inspection,  similar  difficulties  are  evident  in  the  flux  and  source  moment 
expressions  presented  in  the  previous  section.  The  conversion  of  each  equation  to  a 
compact,  computable  form  is  necessary  to  insure  that  any  iteration  procedure  is  sound. 
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i.  Walters  Functions 


Walters  introduced  a  recursive  set  of  exponential  functions  that  implicitly  occur  often 
in  solutions  to  the  Boltzmann  transport  equation  (Walters,  1981:115  and  1986:192-196). 
The  first  three  of  these  expressions,  used  by  Mathews  with  a  lower  case  "p"  (to  avoid 
confusion  with  Legendre  polynomials),  with  z  >  0  are 


,  x  l-exp(-z) 

P°(z)=  z 

(53) 

,  s  l-Po(z) 

P  |(z)  = 

(54) 

l-2p,(z) 

P2(z)~  z 

(55) 

Terms  where  n  >  1  can  be  summarized  by  the  general  Walters  function  forward  iteration 
formula 


\-npn_x{z) 
Pn  (*)  = - ; - 


(56) 


Taking  the  limit  as  the  argument  of  equation  (56)  approaches  zero  yields  a  general 
expression  to  be  used  when  z  =  0: 


pM 


i 

n  +  1 


(57) 


Equations  (53),  (54),  and  (55)  are  plotted  in  Figure  4,  along  with  1/z  and  exp(-z),  with 
z  >  0,  for  comparison. 
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Figure  4.  Functions  1/z,  exp(-z),  and  Walters  Functions  p0(z),  Pi(z),  and  p2(z) 


As  described  by  Mathews,  when  z>n,  equation  (56)  can  accurately  be  used  to 

successively  compute  the  next  higher  term,  up  to  the  level  required  in  a  computation. 
Again,  this  is  known  as  forward  iteration.  If  z  <  n,  the  forward  iteration  process  becomes 
increasingly  less  accurate,  losing  precision  with  each  iteration.  To  overcome  this, 
Mathews  suggests  that  backward  iteration  should  be  performed  by  rearranging  equation 
(56),  so  that 


Pn- .00  = 


1  -zpn(z) 
n 


(58) 


which  is  the  form  for  backward  iteration.  Remarkably,  if  n  is  commenced  at  a  high 
enough  value,  any  initial  guess  in  [0,1]  can  be  used  to  begin  the  backward  iteration 
procedure,  whereupon  an  exact  solution  for  the  Walters  function  can  result  The  worse 
the  initial  guess,  the  higher  the  order  n  must  be  (Mathews,  1990:431-432).  The  error  of 
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any  initial  guess  p  can  be  no  greater  than  the  value  of  the  initial  guess  if  we  use 
1/(2 (n0  +  1))  as  the  initial  guess  (since  1  l(n0  +  1)  >  /?„(£)  >  0  for  £  >  0).  For  a  Walters 
backward  iteration  scheme  that  begins  at  n0,  the  absolute  error  in  the  nth  Walters  function 
using  backward  iteration  is 


AbsErrpn(z) 


n0\ 


and  the  relative  error  is 


(59) 


Errp„(z) 


AbsErrpn(z) 
Exact  pn(z) 


(60) 


A  relative  error  of  less  than  2.0  x  10  17  can  be  realized  for  n  =  0to3  in  backward  iteration 


if  Pr(:)=  1/56  is  used  as  an  initial  guess  for  n0  =  27. 

For  all  transport  calculations  treated  here,  the  highest  order  Walters  functions  required 
for  any  of  the  candidate  methods  in  Table  3  is  p}(z).  If  the  truncated  integer  portion  of  z 
minus  one  is  less  than  or  equal  to  three,  backward  iteration  is  implemented  until  the  order 
of  the  Walters  function  n  is  equal  to  the  truncated  integer  part  of  z.  Forward  iteration  is 
used  thereafter  or  when  the  truncated  integer  part  of  z  minus  one  exceeds  three.  For 
example,  suppose  z  =  2.5.  The  truncated  integer  portion  of  z  less  one  is  1  (the  truncated 
integer  portion  of  2.5  is  2).  Since  1  <  3,  then  backward  iteration  is  used  to  compute 
py{ 2.5)  and  p2( 2.5),  and  forward  iteration  is  used  to  compute  p0( 2.5)  and  p ,(2.5). 
("Walters"  functions  are  used  for  all  the  methods  except  DD  and  LD.  The  linear 
discontinuous  scheme  is  algebraically  equivalent  to  the  linear  nodal  (and  hence  linear 
characteristic)  scheme  when  a  Pade’  (1,2)  approximation  exp(-x)  *  (6  -2x)/(6+x(x  +  4)) 
is  used  in  linear  nodal;  this  Pade’  approximation  is  implemented  when  LD  is  executed, 
and  therefore  use  of  Walters  functions  is  not  necessary). 


38 


If  negative  arguments  are  placed  into  the  Walters  functions,  values  begin  to  increase 
exponentially,  and  vital  digits  can  be  lost  in  successive  calculations  of  large  positive 
numbers.  To  avoid  this,  the  original  Waiters  functions  are  easily  transformed  to  better 
accommodate  negative  arguments  using  combinations  of  the  original  Walters  functions  to 
any  order  desired.  The  transformations  for  equations  (53),  (54),  and  (55)  for  n  =  0to2  are 
(for  i  >  0): 


A)(-z)  =  [Po(z)]exP(z) 

(61) 

px{-z)  =  [p0(z)  -/>,(2)]exp(z) 

(62) 

p2(-z )  =  [p0(z )  -  2px  (z )  +  p2(z )]  exp  (z ) 

(63) 

Using  positive  arguments  and  the  proper  iteration  approach  (forward  or  backward),  stable 
values  can  be  achieved.  Care  must  be  exercised  in  using  these  formulations;  subtraction 
of  nearly  equal  numbers  can  occur  in  some  instances,  resulting  in  a  loss  of  numerical 
precision.  Formulation  of  (61),  (62),  and  (63)  in  truncated  Maclaurin  expansions  of 
suitable  order  is  recommended  in  such  cases.  These  expansions  were  not  required  for 
problems  tested  here. 

ii.  Flux  Equations 

Implementing  the  Walters  functions  into  equations  (44)  to  (48)  results  in  a  very 
compact  set  of  equations.  The  right  boundary  angular  flux  in  equation  (44)  is  equivalent 
to 


V*=Vz.exp[-e]  + 


a  Ax 


p0(P-e)exp[-e] 


(64) 
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where  £  is  the  optical  thickness  for  a  given  p  as  in  equation  (12),  and  p=(-&Ax).  The 

coefficients  a  and  b  are  from  the  assumed  source  function  (equation  (32)).  Similarly,  the 
zeroth  and  first  moments  of  the  angular  flux  are 


and 


Va=VJ>o(p)  + 


a  Av 
P£ 


[P0(P)-Po(P-e)exp[-£]] 


(65) 


Vx=3i|/Le[p2(e)-Pi(e)3+  3aAlP  [p2(p) -/>,$)] 

"2  A 

+-^~r  [2Po(P)-[e  +  2]p0(3-e)exp[-£]] 
p£ 


(66) 


Since  negative  and  positive  argument  forms  for  Walters  functions  are  available  using 
both  forward  and  backward  iteration  schemes  (as  applicable),  no  specific  treatment  of 
terms  using  (P  -  £)  are  typically  required.  Also,  although  z(p2(z)-p,(z))  is  equivalent  to 
(/?„(:)-  2 p,(z)),  the  first  formulation  is  used  in  equation  (66)  (and  wherever  applicable). 
This  is  because  as  ;  — >  0,  (p2(z)  ~Pi(:))  approaches  -1/6,  while  the  other  term  is 
ill-conditioned  and  goes  to  zero.  Further,  although  division  by  £  seems  undesirable, 
these  are  the  most  practical  forms  for  these  equations.  In  the  event  that  optical 
thicknesses  become  extremely  thin  to  the  point  of  causing  catastrophic  cancellations  or 
overflows,  Maciaurin  expansions  of  terms  involving  £  could  be  implemented.  However, 
numerical  testing  of  equations  (64),  (65),  and  (66)  as  presented  with  0.003  revealed  no 
difficulties.  Moreover,  since  the  exponential  characteristic  scheme  is  designed  to  achieve 
high  accuracy  using  coarse  cells,  use  of  an  extremely  fine  mesh  is  unnecessary. 
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iii.  Source  Equations  and  Root  Solvine 


In  order  for  the  exponential  characteristic  method  to  yield  a  solution  to  the  Boltzmann 
transport  equation,  very  accurate  values  for  the  coefficients  a  and  b  are  required.  As 
mentioned  earlier,  equations  (47)  and  (48)  relate  SA  and  Sx  to  the  source  functions 
(equations  (18)  and  (19)).  Since  SA  and  Sx  are  known,  root  solving  can  be  performed 
using  equations  (47)  and  (48)  to  determine  a  and  b  for  each  spatial  cell  in  a  problem. 
Again  making  use  of  the  Walters  functions,  the  following  equations  result  from  a 
manipulation  of  equations  (47)  and  (48): 


SA=ap0(V>) 

(67) 

S,=  [3a(3(p2(p)-p,(p))] 

(68) 

Equation  (67)  can  be  solved  directly  for  a  to  give 

SA 

d  Po((3) 

(69) 

Using  this  and  further  manipulating  equation  (68)  results 

in  an  equation  that  requires  root 

solving  for  (3,  which  indirectly  yields  b ,  since  (3  =  -b Ax  : 

Pi(P) 

Po(P)  ''° 

(70) 

where 


(1-P,) 
r°  2 

and 


(71) 
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(72) 


_SX 
Px  3Sa 

Note  that  no  absolute  values  are  specified  for  equation  (72),  which  differs  from  equation 
(29)  in  the  adaptive  methods  (Mathews,  1990:419-457).  Equation  (70)  can  be  expressed 
in  an  expanded  form  using  exponentials: 


i4+ 


1 


P  exp[p]-l 


-rn 


=0 


A  plot  of  the  function  /?,((3)//?0(P)  is  given  in  Figure  5. 


(73) 


10  -5  0  5  10 

B 

Figure  5.  Ratio  of  p,((3)  to  p0((3)  Walters  Functions 

An  efficient  root  solving  scheme  is  required  to  make  the  exponential  characteristic 
scheme  competitive  with  other  discrete  ordinates  methods.  Newton’s  method  is  preferred 
to  determine  P  for  each  cell,  as  it  has  a  second  order  of  convergence.  One  difficulty  in 
using  Newton's  method  is  that  it  requires  a  reasonable  initial  guess  for  the  root  to 
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guarantee  convergence.  Upon  inspection  of  Figure  5,  note  that  /?,(P)//?0(P)  is  comparable 
to  an  arctan(P)  function.  Normalizing  the  arctangent  function  so  that  it  carries  the  proper 
limits  results  in 


P»(P) 

Po(P) 


^+~arctan[g([})] 
l  71 


(74) 


where  g((3)  is  some  function  which  adjusts  the  curvature  of  the  equation  to  yield  an  exact 

result.  Solving  for  g((3)  in  equation  (74)  and  expanding  the  result  in  a  Maclaurin  series  to 
first  order  yields 


*(P)- 


(#) 
1 12  J 


(75) 


Placing  this  approximation  for  g(P)  back  into  equation  (74)  yields  a  good  approximation 

for  P,(P)/p0(P): 


P  i(P) 
Po(P) 


1  1 

-  + -arc  tan 

2  Jt 


12 


(76) 


J  J 


Replacing  p,(p)/po(P)  in  equation  (70)  with  equation  (76)  and  solving  the  result  for  (3 


yields: 


P« 


(77) 


Note  that  adding  higher  order  terms  in  the  Maclaurin  expansion  for  g([3)  forces  one  to 

solve  for  a  root  of  a  quadratic,  which  defeats  the  purpose  of  finding  a  simple  first  guess 
for  p.  Equation  (77)  provides  an  excellent  first  guess  to  the  root  P  when  0.23  <  r0  <  0.77, 
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with  a  worst  case  absolute  error  of  «  0.3  at  either  endpoint,  becoming  nearly  exact  as 
r0  — »  where  p  — >  0.  A  hyperbolic  tangent  can  also  be  constructed  to  approximate  a  first 
guess  for  the  root;  however,  the  formulation  using  a  tangent  function  is  more  accurate. 

The  tangent  approximation  (equation  (77))  for  p  can  be  used  for  the  entire  range  of  r0, 

which  asymptotically  approaches  one  for  positive  P  and  zero  for  negative  p.  Thu.,,  p 
extends  to  infinity  in  both  the  positive  and  negative  directions.  A  first  guess  better  than 
that  from  equation  (77)  when  r0  >  0.77  is  obtained  by  observing  the  limits  of  the  actual 
equation  (73)  for  large  positive  values  of  p.  When  r0  >  0.77,  then 

Similarly,  when  rQ  <  0.23, 


ft  1 

p  =  --  (79) 

ro 

These  approximations  are  also  very  good,  providing  a  root  to  within  ~  0.3  where  they 

take  over  from  equation  (77),  becoming  nearly  exact  as  r0  approaches  zero  or  one, 
respectively. 

Using  equations  (77),  (78),  or  (79)  as  appropriate  for  an  initial  guess  of  p,  the 

following  equations  for  Newton’s  method  yield  roots  within  5  iterations  at  an  absolute 
tolerance  of  10“15  (determined  by  numerical  testing). 


G(P„,r0) 

G'(P„,r0) 


(80) 


where 
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G(Pn,r0)  = 


1_1  1 
P*  +  (expfpj-  1) 


-rn 


and 


(81) 


G'(P/o)= 


ft  (exp[pj  - 1)  (exp[p„]-l)2 


(82) 


As  of  this  writing,  the  success  of  the  exponential  characteristic  scheme  in  slab 
geometry  (as  is  shown  in  later  sections)  has  prompted  Minor  to  begin  development  of  the 
EC  method  in  two  dimensional  rectilinear  geometry.  Since  the  task  of  root  solving  in 
each  cell  constitutes  the  bulk  of  computational  effort  in  the  EC  scheme,  any  substantial 
improvement  in  root  solving  can  significantly  reduce  the  execution  time  required  for  EC. 
Minor  has  investigated  alternative  root  solving  schemes  to  improve  computational 
efficiency.  Using  equation  (71)  and  solving  equation  (70)  for  p,  and  expanding  the 
numerator  and  denominator  of  the  result  in  a  Maclaurin  series.  Minor  found  (in 
preliminary  studies)  that  the  number  of  iterations  required  for  root  solving  might  be 
reduced  by  at  least  one  iteration.  This  requires  use  of  the  Maclaurin  expansion  over 
defined  sub-domains  wherein  varying  numbers  of  higher  order  terms  are  carried  to  obtain 
double  precision  accuracy.  Minor  is  also  investigating  the  utility  of  a  table-interpolation 
of  Figure  5  to  obtain  a  first  guess  for  root  solving  (Minor,  1991). 


F.  EC  Source  Function  Behavior 

In  order  to  cast  the  distribution  functions  and  moments  into  Walters  functions,  the 
quantity  P  =  -b  Ax  was  introduced,  where  b  is  constant  from  the  source  function 
(equation  (32)).  When  P  is  negative,  b  >  0  and  the  first  source  moment  Sx  is  positive. 
Conversely,  a  positive  value  of  P  indicates  5,  is  negative.  It  is  easily  seen  that  in  a 
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symmetrical  problem,  both  positive  and  negative  values  of  (5  are  encountered, 
emphasizing  the  need  to  treat  Walters  functions  using  negative  arguments  (see  Section 
IV.E.i.). 

Mathematically,  the  exponential  characteristic  scheme  is  robust  due  to  the  unique 
behavior  of  the  source  function  s(x)  =  a[exp(fcr)].  This  function  appears  to  be  a  natural 
choice  if  faced  with  solving  deep  penetration  transport  problems.  If  b  is  identically  zero, 
then  all  equations  describing  flux  distributions  and  moments  default  to  the  step 
characteristic  method.  This  is  useful  in  problems  that  result  in  absolutely  flat  flux 
profiles.  In  the  case  where  ±b  is  very  small  but  non-zero,  as  in  a  region  with  a  virtually 
flat  flux  profile,  the  source  function  can  be  approximated  by  the  truncated  Maclaurin 
expansion 


5(.v)=a(l  +bx)  (83) 

Using  this  approximation  in  place  of  equation  (32)  to  derive  the  flux  distributions  and 
moments  in  equations  (42)  to  (48),  the  exponential  characteristic  method  simplifies  into 
an  algebraic  equivalent  of  the  linear  characteristic  method. 
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V.  Program  Development  and  Validation 

Program  development,  as  discussed  here,  consisted  of  three  phases:  standardization, 
optimization,  and  computer  implementation  of  the  exponential  characteristic  spatial 
quadrature  method.  Standardization  dealt  with  adapting  existing  codes  using  candidate 
comparison  methods  (DD,  DDF,  LD,  LC,  SA,  and  LA)  into  one  complete  code,  so  that 
all  quantities  could  be  computed  using  standard  variables  and  methods.  To  optimize  the 
single  multi-method  code,  several  steps  were  taken,  including  the  installation  of  timing 
sequences,  initializations,  and  moment  and  overall  balance  checks,  as  well  as  the 
streamlining  of  data  reporting  and  storage.  Finally,  the  exponential  characteristic 
method,  using  computable  forms  for  equations  derived  in  the  previous  section,  was 
implemented  into  the  multi-method  code.  As  for  program  validation,  all  methods  were 
continuously  validated  and  checked  for  accuracy  following  each  major  code 
modification. 

A.  Standardization  and  Optimization 

Several  codes  were  supplied  by  LCDR  K.  Mathews,  Ph.D.  for  solving  slab  geometry 
discrete  ordinates  problems.  Separate  codes  using  diamond  difference  (with  a  fixup 
option),  linear  discontinuous,  step  adaptive,  and  linear  adaptive  spatial  quadratures  were 
supplied.  A  linear  characteristic  routine  was  constructed  by  modifying  the  linear 
adaptive  model.  Because  each  code  was  slightly  different  in  the  way  values  were 
computed,  every  effort  was  made  to  standardize  all  calculations  into  a  single, 
multi-method  code.  A  standardized  code  was  successfully  implemented  and  allows  the 
computational  effort  for  any  spatial  quadrature  scheme  to  be  determined  directly  (for  a 
specific  problem)  from  the  execution  time  required  to  achieve  a  converged  solution.  All 
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execution  times  were  computed  based  on  real  clock  inquiries  from  the  computer 
operating  system.  (To  avoid  single  system  multi-tasking  timing  conflicts,  no  transport 
solutions  were  obtained  under  the  Microsoft  Windows  3.0  PC  environment). 

The  boundary  conditions  treatable  in  this  multi-method  code  include  a  choice  of  a 
vacuum,  a  symmetric  albedo,  or  a  grey  albedo  at  each  slab  face.  Types  of  incident 
currents  available,  also  at  each  slab  face,  are  a  Lambertian,  an  isotropic  surface  source,  or 
a  coliimated  beam  (at  any  angle  of  incidence  to  the  slab).  These  boundary  conditions  and 
incident  currents  are  discussed  further  in  Appendix  B.  All  iterations  begin  at  the  left  slab 
edge  with  positive  angles  and  progress  from  left  to  right,  then  from  right  to  left  with 
negative  angles.  Fluxes  and  currents  are  all  initialized  to  zero  as  each  different  spatial 
scheme  is  applied  to  a  problem.  The  balance  equations  (equations  (25)  and  (26))  are 
verified  for  each  cell  during  each  iteration.  The  left  and  right  sides  of  these  equations  are 
compared,  and  the  resulting  maximum  differences  are  reported.  Upon  convergence  to  a 
specified  relative  tolerance  (1CT5),  angular  quadratures  are  performed  and  scalar  fluxes 
and  currents  are  globally  computed.  Region  values  are  then  calculated  by  folding 
together  cell  values.  The  overall  balance  equation  (equation  (33))  is  verified  in  a  manner 
similar  to  that  for  equations  (25)  and  (26). 

All  input  files  for  any  problem  can  be  set  up  in  a  standard  format  specified  in  the 
MAKEIN1D  input  code.  Due  to  unique  differences  involved  in  treating  the  source 
function  using  fixed  point  and  successive  scatter  methods,  separate  multi-method  codes 
are  used  for  each  technique.  The  name  assigned  to  the  fixed-point  multi-method  code  is 
SN1D-ALL,  while  the  name  assigned  to  the  successive  scatters  version  is  SS ID-ALL. 
Full  featured  output  options  are  available,  including  ASCII  data  files  which  facilitate 
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plotting  of  solutions.  In  addition,  a  separate  error  and  order  of  convergence  processing 
algorithm,  PR0CSN02,  was  developed  to  deal  with  the  results  from  a  large  collection  of 
transport  solutions  for  any  given  problem. 

B,  Implementation  of  EC  Spatial  Quadrature 

Following  all  of  the  modifications  required  for  a  smooth,  multi-method,  standardized 
code,  only  one  step  was  required  to  implement  exponential  characteristic  spatial 
quadrature.  This  step  consisted  of  developing  the  subroutine  which  computed  the  angular 
fluxes,  based  on  source  moments,  for  each  mesh  cell. 

Naturally,  the  subroutine  for  the  exponential  characteristic  method  had  to  include 
code  for  root  solving  to  determine  the  coefficients  a  and  b  in  equation  (32).  Following 
an  initial  guess  for  p  using  either  equations  (77),  (78),  or  (79),  root  solving  using 
equations  (80),  (81),  and  (82)  was  performed.  Performance  studies  were  conducted  to 
determine  the  optimum  relative  tolerance  used  in  the  Newton  iteration  loop.  A  relative 
root  tolerance  of  10"5  was  selected.  This  value  provided  a  root  that  was  accurate  enough 
for  problem  convergence  to  10“5,  yet  allowed  roots  to  be  found  rapidly.  In  addition,  flux 
moments  computed  from  roots  calculated  using  this  tolerance  satisfied  the  balance 
equations  (See  equations  (25),  (26),  and  (33))  as  well  as  the  other  methods.  Typically, 
between  two  and  five  iterations  were  required  to  find  a  root.  Decreasing  the  root  solving 
tolerance  below  10~5  causes  EC  to  require  increased  execution  time.  As  shown  in  test 
problems  presented  in  subsequent  sections,  use  of  a  relative  tolerance  of  10~5  offered 
excellent  relative  performance. 

Following  the  determination  of  the  root  P  =  -bAx,  Walters  functions  using  P  and  £ 
(the  angular  optical  thickness)  were  calculated.  If  needed,  the  negative  arguments  for  the 
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Walters  functions  were  treated  using  equations  (61),  (62),  and  (63),  and  all  required 
fluxes  and  flux  moments  were  calculated  using  equations  (64),  (65),  and  (66).  No  other 
functions  were  necessary  inside  the  exponential  characteristic  spatial  quadrature  step. 

The  subroutine  which  performs  these  tasks  is  "StepEC"  in  the  source  code  for 
SN1D-ALL  in  Appendix  D. 

C.  Validation 

The  multi-method  slab  geometry  codes  were  verified  by  comparing  solutions  of 
conventional  methods  with  independent  versions  of  the  originally  supplied  codes  for  DD, 
LD,  SA,  and  LA.  An  additional  validation  was  performed  using  one  independent  slab 
geometry  code  employing  the  DD  scheme.  Problems  used  to  verify  all  solutions 
consisted  of  a  selection  of  optically  thin  single  and  multiple  region  problems  (with 
various  scattering  and  absorption  properties).  All  methods  yielded  essentially  identical 
numerical  solutions  to  these  problems.  Solutions  to  these  validation  problems  using  the 
exponential  characteristic  method  agreed  with  those  derived  using  conventional  methods. 


50 


VI.  Testing  and  Evaluation 


While  many  slab  geometry  transport  problems  were  solved  during  the  development  of 
the  exponential  characteristic  scheme,  only  the  results  from  a  few  select  test  cases  are 
presented.  Clearly,  there  are  an  infinite  number  of  scenarios  where  the  performance  of 
the  EC  method  versus  other  spatial  quadrature  methods  might  be  investigated.  However, 
problems  put  forth  here  best  illustrate  the  strengths  and  weaknesses  of  the  exponential 
characteristic  method  discovered  to  date.  Each  test  problem  was  solved  using  both  fixed 
point  and  successive  scatter  iterations  of  the  source  function  with  various  mesh  sizes 
corresponding  to  a  range  of  optical  thicknesses  a  Ax.  Limiting  assumptions,  comparisons 
to  other  methods,  error  analyses,  and  computational  costs  and  efficiencies  are  treated  as 
described  in  Section  III. 

A.  Test  Case  I;  Thick  Absorber 

To  initially  challenge  the  exponential  characteristic  method  and  evaluate  its  utility,  a 
single  region  with  an  optical  thickness  of  sixteen  (used  by  Mathews  in  his  evaluation  of 
the  linear  adaptive  method)  was  used  as  the  first  test  case  (Mathews,  1990:  444-445).  A 
schematic  of  this  problem  is  presented  in  Figure  6. 

The  source  on  the  left  boundary  (an  ideal  surface  of  isotropic  emitters)  has  an 
angular  flux  distribution  inverse  to  |i,  corresponding  to  large  fluxes  at  steep  incident 
angles  (See  Appendix  B).  Augmented  by  a  small  amount  of  scatter,  the  scalar  flux  on  the 
left  boundary  for  the  reference  solution  is  more  than  3.0.  The  region  average  scalar  flux, 
as  expected,  is  reduced  to  nearly  l/45th  of  this  value,  reaching  a  value  of  =  10“*  at  the 
right  boundary.  The  reference  solution  was  computed  with  the  linear  characteristic 
method  using  256  cells  over  the  single  region.  AH  methods  performed  well  at  the  left 
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Vacuum 


Figure  6.  Schematic  for  Test  Case  1 


boundary  where  current  enters  the  slab.  However,  of  real  interest  is  how  accurately  the 
various  methods  (listed  in  Table  3)  computed  the  region  average  and  right  boundary 
(x  =  16)  scalar  fluxes.  Error  norms  for  each  spatial  quadrature  method  for  the  region 
average  scalar  flux  are  presented  in  Figure  7.  This  figure  clearly  shows  that  the 
exponential  characteristic  (EC)  method  outperformed  all  of  the  others,  achieving  results 
for  average  flux  within  less  than  0.2  percent  of  the  reference  solution  using  an  optical 
thickness  of  4  (corresponding  to  only  4  mesh  cells).  Also  note  that  EC  performed  nearly 
as  well  using  only  a  single  sixteen  mean  free  path  cell. 

The  error  norms  for  the  right  boundary  flux  are  presented  in  Figure  8,  where  again, 
the  exponential  characteristic  method  yielded  much  more  accurate  results  using  thick 
cells.  Observe  in  Figure  8  that  using  an  optical  thickness  cAx  =  4,  the  error  of  the  EC 
right  boundary  flux  was  less  than  2  percent,  while  for  LA,  the  nearest  competitor,  the 
error  was  greater  than  10  percent.  On  the  coarsest  mesh,  EC  performed  slightly  better 
than  LA.  These  results  are  typical  of  those  found  for  currents  and  first  flux  moments  as 
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Test  Case  i  Optical  Thickness  per  Cell 

Fixed  Point  Source  Iteration 


Figure  7.  Test  Case  1  Error  in  Region  Average  Scalar  Flux 


well.  Note  that  while  the  diamond  difference  (DD)  scheme  converged  rapidly  as 
oAt  — »  1  for  the  region  average  flux  in  Figure  7,  it  was  severely  in  error  at  this  same 
optical  thickness  at  the  right  boundary,  as  shown  in  Figure  8. 

In  terms  of  average  pointwise  error  Q,  computed  using  equation  (37),  the  EC  scheme 
was  almost  four  times  more  accurate  than  its  nearest  competitor,  LA.  Also,  as  can  be 
discerned  from  Figures  7  and  8,  the  EC  method  maintained  a  bounded  error  using 
optically  thick  cells  in  this  problem.  As  the  mesh  became  highly  refined,  EC  showed  at 
ica.M  4th  order  convergence.  (Convergence  appeared  to  be  5th  order  in  the  limit  of  very 
thin  cells;  this  was  not  confirmed  due  to  limitations  in  root  solving,  which  required  a  shift 
to  LC  when  using  extremely  fine  cells,  as  discussed  in  the  next  section). 

Mesh  size  ratio  (MSR),  discussed  in  Section  m.C,  is  useful  in  illustrating  a 
normalized  mesh  size  relative  to  the  number  of  cells  used  by  the  exponential 
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Figure  8.  Test  Case  1  Error  in  Right  Boundary  Scalar  Flux 


characteristic  scheme  to  achieve  an  average  pointwise  error  Q.  Figure  9  is  a  graph 
depicting  the  relative  MSR  for  this  single  thick  absorber  problem.  (MSR  values  were 
computed  using  linear  interpolation  of  average  pointwise  error  results  over  optical 
thicknesses  typical  of  those  shown  in  Figures  7  and  8).  Observe  in  Figure  9  that  when  4 
mesh  cells  were  used  for  EC  quadrature  (yielding  an  error  Q  less  than  1.5  percent),  the 
MSR  for  linear  adaptive  quadrature  was  2.5.  This  means  that  for  LA  to  achieve  a  average 
pointw  ise  error  less  than  1.5  percent,  two  and  one  half  times  the  number  of  cells  used  by 
EC  were  necessary  (10  cells). 

Figure  9  further  demonstrates  that  significantly  more  cells  are  required  by  the  other 
methods  for  an  equivalent  average  error,  DD  and  DDF  would  require  such  an  abundance 
of  cells  to  obtain  a  calculation  equivalent  to  EC  that  they  offered  no  competition.  Note 
from  the  figure  that  both  LD  and  SA  appeared  to  level  off  with  MSRs  of  nearly  5  and  3, 
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Fixed  Point  Source  Iteration 


Cells  Used  for  EC 


Figure  9.  Test  Case  1  Mesh  Size  Ratio 

respectively,  while  LC  and  LA  merged  to  an  MSR  value  of  1.7  as  the  mesh  was  refined. 
Obviously,  for  this  single  thick  absorber,  EC  is  superior.  In  comparing  execution  times 
required  to  reach  an  equivalent  error  Q ,  the  EC  method  offered  an  advantage  when  up  to 
five  cells  were  used  (with  an  error  Q  of  one  percent).  However,  exec  mon  times  for 
single  region  absorption  problems  were  so  low  that  actual  comparisons  are  made  only  in 
later  problems  with  at  least  two  regions. 


Solutions  for  this  problem  using  successive  scatters  of  the  source  function  yielded 
similar  results.  Overall,  errors  computed  for  successive  scatters  were  10  percent  lower 
for  exponential  characteristic  quadrature,  and  remained  mostly  unchanged  for  all  of  the 
other  methods.  Apparently,  EC  benefits  from  tracking  each  scattered  flight  of  neutrons  in 
this  thick  absorber,  possibly  because  it  more  accurately  represents  the  source  moment 
behavior  of  each  flight  with  the  exponential  source  function.  In  any  event,  exponential 
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characteristic  has  proven  to  be  a  very  robust  scheme  in  test  case  1. 


B.  Test  Case  2:  Moderate  Scatter  and  Source 

Because  the  exponential  characteristic  method  showed  much  promise  in  the  first  one 
region  absorber  problem,  two  different  regions  were  specified  for  the  second  test  case.  A 
schematic  of  this  problem  is  given  in  Figure  10.  The  first  region  is  a  sixteen  mean  free 
path  thick  absorber  with  moderate  scattering  and  is  coupled  to  a  second  equally  thick 
absorbing  region  containing  a  source.  The  left  boundary  is  a  vacuum,  and  the  right 
boundary  is  an  isotropic  surface  source  with  a  vacuum. 


Vacuum 


Isotropic 


Optical  Thickness  =  16.0  per  Region 


Figure  10.  Schematic  for  Test  Case  2 


With  moderate  scattering  in  the  left  region,  a  modest  source  in  the  right  region,  and 
an  isotropic  surface  source  on  the  right  boundary,  the  flux  profile  in  the  second  (source) 
region  should  be  almost  flat.  If  this  is  the  case,  the  constant  b  in  the  EC  source  function 
equation  (32)  will  be  very  small,  and  the  EC  method  can  be  expected  to  perform  like  the 
LC  method,  as  discussed  in  Section  IV. F. 
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Test  Case  2 

Fixed  Point  Source  Iteration 

Optical  Thickness  of  8,  2  Cells  per  Region 


Distance,  cm 


Figure  11.  Test  Case  2  Cell  Boundary  Scalar  Flux 


The  cell  bo^.idary  scalar  fluxes  for  this  problem  are  shown  in  Figure  1 1  using  two 
cells  per  region,  equivalent  to  an  optical  thickness  aA.r  =  8  in  each  cell.  Results  of  region 
average  scalar  fluxes  are  reported  in  Table  4  for  this  same  optical  thickness.  As  in  the 
first  problem,  a  reference  solution  was  computed  with  linear  characteristic  spatial 
quadrature  with  256  cells  per  region.  Observe  that  the  flux  was  reasonably  flat  in  the 
right  region  between  x  =  16  and  x  =  32,  as  anticipated. 

Moving  to  the  left  from  the  right  edge  at  x  =  32  to  the  cell  interface  in  the  middle  of 

region  2  at  .v  -  24.  all  methods  appeared  to  yield  effectively  the  same  result  except  for  the 
diamond  difference  method,  which  plunges  to  a  meaningless  negative  result.  Between 
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Table  4 

Test  Case  2  Region  Average  Scalar  Fluxes  for  o  At  =  8 


Method 

Region  1 

Region  2  <|>A 

DD 

2.8627E-02 

5.7524E-01 

DDF 

2.835  IE-02 

6.0358E-01 

LD 

7.9753E-02 

6.0516E-01 

LC 

1.6436E-02 

6.0545E-01 

SA 

1.4972E-02 

6.0653E-01 

LA 

1.5143E-02 

6.0617E-01 

EC 

1.4876E-02 

6.063  IE-01 

Ref 

1.4646E-02 

6.0550E-01 

x  =  16  and  x  =  24,  DDF  and  especially  DD  were  noticeably  in  error.  At  x  =  16  where 
regions  1  and  2  meet,  small  differences  in  flux  were  reported  by  the  various  methods,  as 
shown  in  the  exploded  portion  of  the  graph  (note  how  closely  EC  tracked  along  the 
reference  solution).  In  the  left  region  between  x  =  0  and  x  -  16,  a  wide  variety  of 
solutions  resulted,  indicating  severe  numerical  diffusion  from  the  use  of  optically  thick 
cells.  Remarkably,  the  flux  solution  for  the  exponential  characteristic  scheme  agreed 
almost  precisely  with  the  reference  solution  throughout  both  regions.  As  in  test  case  1, 
LA  offered  the  next  best  solution  to  EC,  followed  closely  by  SA. 

Results  of  average  pointwise  and  maximum  observed  errors  are  listed  in  Table  5  for 
an  optical  thickness  of  aAr  =  8.  At  this  optical  thickness,  the  EC  method  yielded  an 
average  pointwise  error  Q  less  than  0.04,  while  this  error  for  LA  was  in  excess  of  0.26. 
For  the  exponential  characteristic  scheme,  the  maximum  error  observed  at  any  boundary 
or  region  value  for  oAt  =  8  (including  all  fluxes,  currents,  and  moments)  for  this  problem 
occurred  in  the  solution  for  the  net  current  at  the  left  boundary  (x  =  0),  which  was  0.14. 
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The  maximum  error  for  the  LA  scheme  was  found  in  the  solution  for  the  flux  on  the  left 
boundary,  which  was  1.36.  This  further  demonstrates  the  ability  of  the  EC  scheme  to 
yield  a  reliable  solution  over  a  very  coarse  mesh. 


Table  5 

Test  Case  2  Error  Results  for  ctAx  =  8 


Method 

Average  Error 

Q 

Maximum 
Observed  Error 

Location  of  Maximum 
Observed  Error 

DD 

0.9429 

2.0000 

Left  Bdy  <j> 

DDF 

0.7580 

2.0000 

Left  Bdy  0 

LD 

0.6742 

2.0000 

Region  1  0, 

LC 

0.4597 

1.9823 

Left  Bdy  0 

SA 

0.2853 

1.4231 

Left  Bdy  0 

LA 

0.2642 

1.3651 

Left  Bdy  0 

EC 

0.0388 

0.1409 

Left  Bdy 

A  difficulty  arose  in  the  EC  method  for  test  case  2  when  optical  thicknesses  dropped 
below  unity  during  calculations  in  the  right  source  region.  Analysis  of  intermediate 
iteration  values  revealed  that  the  first  moment  of  the  source  was  very  small  when 
c\x  <  ~1.  The  ratio  of  the  first  source  moment  Sx  to  the  zeroth  moment  SA  was  less  than 
6.0  x  KT6  for  some  cells.  The  value  of  [3  =  -frAx  was  almost  zero,  and  the  number  of 
root  solving  iterations  required  by  Newton’s  method  (See  Section  IV.E.iii)  became 
unfavorably  large,  even  for  modest  comparison  tolerances.  Although  in  these  cases  the 
exponential  characteristic  scheme  is  algebraically  equivalent  to  the  linear  characteristic 
method,  root  solving  is  still  required  to  obtain  (3  and  therefore  b.  When  these 
complications  arose,  the  execution  time  for  the  EC  scheme  became  extremely  large, 
burdened  by  the  root  solving  overhead  progressing  into  thousands  of  iterations  in  some 
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cells,  as  opposed  to  the  typical  two  to  five  normally  required  following  the  initial  guess. 
When  some  very  optically  thin  cells  were  specified,  the  flux  profile  changed  so  little 
across  some  cells  that  no  roots  were  found;  the  Newton  iteration  scheme  oscillated  the 
root  around  zero. 

Since  the  EC  method  only  encounters  difficulty  in  root  solving  when  Sx  is  nearly 

zero  (compared  to  SA ),  and  since  it  is  asymptotically  equivalent  to  the  linear  characteristic 
method  in  these  circumstances  (See  Section  IV. F),  a  switch  to  LC  was  implemented  into 
the  EC  algorithm  when  the  initial  guess  for  P  fell  to  less  than  or  equal  to  1.2  x  10~5.  This 
value  was  found  to  be  a  lower  limit  to  maximize  the  use  of  the  EC  method  and  still  avoid 
root  solving  problems.  The  upper  limit  to  this  trip  point  is  defined  to  be  when  the  LC 
scheme  can  be  actuated  with  no  fixup.  This  limit  is  satisfied  when  P  is  as  large  as  two. 
However,  setting  the  trip  point  to  use  LC  when  P  <  2  sacrifices  the  advantage  of  the  EC 
method  in  many  cells,  and  could  result  in  oscillating  or  failed  convergence.  This  is 
because  the  exponential  source  moments  in  EC  can  be  quite  unlike  the  linear  moments  in 
LC  for  values  of  p  not  close  to  zero,  potentially  resulting  in  iterations  that  flip-flop 
between  EC  and  LC  and  thus  never  converge.  However,  this  difficulty  was  not  observed 
in  test  case  2  when  the  trip  point  for  P  was  set  to  2. 

In  an  effort  to  deteimine  an  optimum  set  point,  limned  tests  using  test  case  2  were 
conducted  with  several  optical  thicknesses,  computed  over  a  range  of  increasing  initial 
guess  setpoints  for  P  to  change  to  LC  (between  1.2  x  10"5  and  2.0).  As  the  set  point  was 
increased,  resulting  solutions  became  less  accurate.  However,  between  P  <  1.2  x  10~5  and 
up  to  P  <  5.0  x  10~\  no  changes  in  any  solutions  were  detected  (using  the  standard 
relative  tolerance  of  10_s),  but  first  moment  balances  (using  equation  (26))  were 
conserved  to  a  much  higher  level  using  P  <  5.0  x  lO^4.  With  a  trip  point  p  <  1 .2  x  10"\  the 
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first  moment  relative  tolerance  was  ~  while  using  p  <  5.0  x  10-4  yielded  a  tolerance 
of  =  10~10.  An  optimum  set  point  of  P  <  5.0  x  10~*  corresponds  to  an  SX/SA  of  2.5  x  1CT\ 
Further,  this  prevents  extraordinarily  high  numbers  of  iterations  for  root  solving  in  EC, 
maximizes  moment  balance  tolerances,  and  actuates  a  switch  to  LC  for  a  few  cases  when 
root  solving  becomes  uniquely  difficult  (for  very  flat  flux  profiles  across  a  cell,  when 
Sx  — »  0).  Under  these  constraints,  LC  never  requires  a  fixup  and  conserves  all  moments. 
Also,  use  of  the  LC  method  for  an  initial  guess  for  P  <  5.0  x  KT*  is  reasonable,  as  the 
exponential  characteristic  scheme  is  asymptotically  equivalent  to  linear  characteristic  in 
this  case.  As  a  result,  solutions  for  test  case  2  were  computed  using  a  setpoint  of 
p<  5.0x1  O'4  switch  to  LC. 

A  plot  of  the  average  pointwise  error  using  a  fixed  point  source  iteration  for  this 
problem  is  given  in  Figure  12.  As  shown,  errors  for  the  EC  method  were  far  lower  than 
for  any  other  scheme,  especially  for  a  coarse  spatial  mesh. 


Figure  12.  Test  Case  2  Average  Pointwise  Error,  Fixed  Point  Source  Iteration 
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Analysis  of  successive  scatter  solutions  to  test  case  2  demonstrated  that  the  switch  to 
the  LC  method  using  P  <  5.0  x  1(T*  was  not  necessary  until  optical  thicknesses  fell  below 
aAv  =  2.0,  compared  to  oAx  =  4.0  when  fixed  point  source  iteration  was  used.  This 
result  is  observed  because  in  successive  scatters,  the  external  source  is  only  present  in  the 
source  function  during  the  initial  flight  of  neutrons  in  the  source  region.  After  the  first 
flight,  the  magnitude  of  the  source  function  comes  from  only  scattering  of  the  first  flight. 
The  source  term  is  not  dominated  by  the  constant  external  source  (as  in  the  fixed  point 
case)  in  successive  flights,  resulting  in  ratios  of  first  to  zeroth  source  moments  that  are 
not  as  close  to  zero.  When  this  occurs,  P  is  not  exceedingly  small.  The  need  to  use  LC 
only  arises  when  the  optical  thickness  is  refined  to  a  degree  that  the  scattered  flux  for  a 
given  flight  of  neutrons  is  very  flat  across  a  cell,  yielding  values  of  SX/SA  and  p  that  are 
small. 


Figure  13.  Test  Case  2  Average  Pointwise  Error,  Successive  Scatters  Source  Iteration 
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The  average  pointwise  errors  Q  for  test  case  2  using  successive  scatters  are  shown  in 
Figure  13.  The  greatest  contribution  to  the  average  error  occurred  for  both  fixed  point 
and  successive  scatters  at*  =  0.  Using  very  coarse  meshes,  the  error  using  successive 
scatters  was  almost  twice  that  of  fixed  point  source  iteration.  This  was  not  observed  in 
the  first  problem  which  had  a  strongly  absorbing  region.  Apparently,  when  very  coarse 
meshes  are  specified  with  more  than  minimal  scattering,  successive  scatters  performs 
worse  than  fixed  point  iteration  of  the  source  function.  At  an  optical  thickness  of  two,  the 
errors  by  either  method  are  nearly  equal.  Below  oAr  =  2.0,  a  benefit  of  using  successive 
scatters  was  realized  with  decreasing  optical  thickness  because  the  source  function  was 
not  dominated  by  the  external  source  term;  using  optically  thin  cells,  average  pointwise 
errors  for  successive  scatters  were  only  =60  percent  of  those  for  fixed  point  iteration. 
Nevertheless,  since  the  EC  scheme  delivers  accurate  solutions  using  a  only  a  few  cells, 
these  results  suggest  that  for  this  problem,  the  additional  five  iterations  required  for  a 
successive  scatters  solution  are  not  worth  the  computational  effort. 

Analogous  to  the  convergence  in  test  case  1,  as  the  optical  thickness  oAt  —>  0,  the 

EC  method  went  to  4th  order  in  test  case  2.  The  convergence  of  the  LC  and  LA  schemes 
trailed  that  of  EC  over  coarser  cells,  each  close  to  one  at  first,  but  also  ending  at  4th  order 
as  crAv  — >  0.  As  expected,  the  diamond  schemes  began  with  very  low  convergence  and 
closed  on  second  order  as  the  optical  thickness  decreased.  The  convergence  of  LD  and 
SA  showed  an  improvement  over  DD  and  DDF.  Overall,  the  EC  method  demonstrated 
the  highest  order  of  convergence  in  the  limit  of  thin  cells. 

The  mesh  size  ratio  (MSR)  for  test  case  2  is  presented  in  Figure  14.  These  results 
are  similar  to  those  presented  in  Figure  9  for  test  case  1.  When  EC  used  2  cells,  an  MSR 
of  nearly  4  for  LC  and  LA  resulted,  indicating  an  equivalent  of  8  cells  for  LC  and  LA. 
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When  4  cells  were  used  with  the  exponential  characteristic  scheme,  the  linear  adaptive 
and  linear  characteristic  methods  required  3  times  as  many  (12  cells)  to  obtain  the  same 
error  Q . 


Figure  14.  Test  Case  2  Mesh  Size  Ratio 


Execution  time  ratios  (ETRs),  computed  using  equation  (39),  are  presented  in  Figure 
15.  Here  the  ETR  for  each  spatial  quadrature  scheme  is  plotted  versus  average  pointwise 
error  Q,  and  ETR  values  are  execution  times  normalized  to  the  time  used  by  the  EC 
scheme.  At  an  error  Q  slightly  less  than  4  percent,  the  EC  scheme  used  2  mesh  cells,  as 
indicated,  yielding  a  corresponding  ETR  of  =2.5  for  LC  and  LA.  Therefore,  LC  and  LA 
required  about  2.5  times  the  execution  time  as  that  consumed  by  the  EC  scheme  to  obtain 
less  than  a  4  percent  error  Q .  Only  when  the  average  pointwise  error  fell  below  0.6 
percent  did  the  execution  time  for  EC  equal  that  for  LC  or  LA  (using  more  than  15  cells 
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Figure  15.  Test  Case  2  Execution  Time  Ratio 


each)  when  =7  cells  were  specified  for  EC.  Thus,  if  memory  storage  is  not  a  concern, 
either  LC,  LA,  or  EC  would  provide  an  average  pointwise  error  of  less  than  0.6  percent  in 
close  to  the  same  execution  time. 

C.  Test  Case  &  Absorber  and  Diffusion 

Clearly,  the  exponential  characteristic  spatial  quadrature  scheme  has  demonstrated 
superior  performance  over  all  of  the  conventional  schemes  tested.  In  test  case  2,  the  left 
region  was  a  modest  scatterer,  and  the  EC  scheme  yielded  the  best  overall  solution,  even 
better  than  LA  or  LC  using  optically  thin  cells.  It  is  known  that  the  diamond  difference, 
linear  discontinuous,  and  linear  characteristic  methods  satisfy  the  diffusion  limit;  as  a 
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diffusive  problem  is  refined,  these  transport  methods  produce  solutions  that 
asymptotically  approach  the  solution  produced  by  the  neutron  diffusion  equation  (Larsen, 
1982:  90-95).  In  diffusive  problems,  Sx  is  typically  much  smaller  than  SA,  and  in  the  limit 
as  the  mesh  is  refined,  the  EC  scheme  is  asymptotically  equivalent  to  the  LC  scheme  as 
Sx  — » 0.  Moreover,  since  the  EC  scheme  performed  better  than  the  DD,  LD,  and  LC 
schemes  in  the  moderately  diffusive  region  of  test  case  2, 1  anticipated  that  the 
exponential  characteristic  scheme  would  satisfy  the  diffusion  limit.  To  test  this 
expectation,  test  case  3  contained  two  regions;  a  strongly  absorbing  region  on  the  left, 
and  a  highly  diffusive  region  with  a  source  on  the  right.  Also,  an  isotropic  surface  source 
was  located  at  the  right  boundary,  and  vacuum  boundaries  were  present  on  either  side.  A 
schematic  of  this  problem  is  given  in  Figure  16. 

Isotropic 


Optical  Thickness  =  16.0  per  Region 
Figure  16.  Schematic  for  Test  Case  3 

The  solutions  to  the  scalar  flux  for  this  problem  are  presented  below  in  Figure  17. 
Only  two  cells  were  used  to  obtain  these  solutions;  the  optical  thickness  was  oAv  =  8.  As 
expected,  the  flux  solutions  for  the  various  methods  were  very  close  as  a  result  of  the 
strong  diffusion  of  neutrons  from  the  right  region.  As  indicated  in  the  figure,  the 


66 


exponential  characteristic  scheme  was  not  as  accurate  in  predicting  fluxes  as  in  other 
problems.  The  exponential  source  function  in  EC  has  more  difficulty  in  modeling  the 
scalar  flux  in  diffusive  regions  than  the  linear  source  function  in  LC.  This  is  expected,  as 
an  exponential  approximation  of  a  concave  downward  flux  profile  is  frequently  more  in 
error  than  for  a  linear  approximation.  Still,  results  for  EC  were  found  to  agree  closely 
with  LC  and  LA. 


Figure  17.  Test  Case  3  Cell  Boundary  Scalar  Flux 
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Curiously,  when  boundary  currents,  moments,  and  scalar  fluxes  are  all  considered, 
the  EC  scheme  yielded  the  lowest  average  pointwise  error  Q.  A  plot  of  this  average  error 
is  given  in  Figure  18.  Note  in  the  figure  that  the  error  Q  for  the  LA  scheme  was  slightly 
greater  than  that  for  EC. 


Tes!  Case  3  Optical  Thickness  per  Cell 

fixed  Point  Source  Iteration 


Figure  18.  Test  Case  3  Average  Pointwise  Error 


Furthermore,  considering  the  maximum  observed  error  in  any  region  or  boundary 
flux,  current,  or  moment,  the  EC  scheme  yielded  very  favorable  results.  A  plot  of  the 
maximum  observed  error  for  problem  3  is  presented  in  Figure  19.  Although  the 
maximum  observed  error  for  LA  and  EC  are  nearly  the  same,  the  average  pointwise  error 
of  the  exponential  characteristic  scheme  was  still  lower  than  that  for  the  linear  adaptive 
method.  Use  of  a  successive  scatters  iteration  of  the  source  function  yielded  solutions 
that  were  slightly  degraded  for  all  methods.  In  most  cases,  the  difference  between  the 
fixed  point  solution  and  the  successive  scatters  solution  was  minor. 
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Test  Case  3  Optical  Thickness  per  Cell 

Fixed  Point  Source  Iteration 

Figure  19.  Test  Case  3  Maximum  Observed  Error 

The  EC  solution  to  this  problem  proved  to  be  in  reasonable  agreement  with  the 
solutions  for  DD,  LD,  and  LC  methods.  Of  these  methods,  EC  yielded  the  lowest 
average  and  maximum  errors  (when  considering  all  computed  flux  and  boundary  values). 
Experience  using  the  exponential  characteristic  method  in  other  extremely  diffusive 
problems  showed  that  EC  performed  better  globally  than  either  DD,  DDF,  or  LD  over 
any  optical  thickness.  In  many  cases,  the  solution  afforded  by  EC  was  only  slightly 
worse  than  that  for  LC  or  LA;  again,  this  is  likely  due  to  the  source  function  in  LC  and 
LA  using  a  linear  approximation  instead  of  an  exponential  as  in  EC.  Sometimes,  EC 
yielded  the  best  overall  solution  by  a  small  margin.  Although  supporting  evidence  is 
extremely  limited,  these  findings  suggest  that  the  EC  method  asymptotically  approaches 
the  thin  cell  diffusion  limit. 
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Analysis  of  the  mesh  size  ratio  for  test  case  3  revealed  that  little  advantage  was 
gained  in  using  the  EC  scheme  over  the  linear  characteristic  or  linear  adaptive  methods, 
as  might  be  expected  after  inspection  of  the  average  pointwise  error  in  Figure  18.  As 
mentioned,  solutions  for  the  EC,  LA,  and  LC  schemes  were  very  close.  When  execution 
time  ratios  were  compared,  the  EC  method  was  at  a  disadvantage,  requiring  more  than 
twice  the  execution  time  of  LA  for  the  same  error.  Thus,  on  the  basis  of  accuracy  and 
execution  time  in  test  case  3,  EC  was  not  the  most  efficient  scheme. 
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VII. 


From  the  results  for  test  case  1,  the  thick  absorber,  the  exponential  characteristic 
scheme  was  superior  to  any  other  method.  This  is  credited  to  the  natural  ability  of  the 
exponential  source  function  to  model  nearly  exponential  attenuation.  Other  spatial 
quadratures  required  at  least  2  to  3  times  the  number  of  cells  that  EC  needed  to  yield  an 
equivalent  solution.  This  reveals  the  potential  savings  in  memory  storage  offered  by  EC 
when  applied  to  deep  penetration  absorbers.  In  two  dimensions,  the  memory  required  by 
EC  may  be  as  little  as  one  tenth  of  that  necessary  for  most  other  methods  for  a  given 
accuracy. 

The  strength  of  the  exponential  characteristic  scheme  was  further  demonstrated  in  test 
case  2.  This  problem  included  a  moderate  scatterer  region  and  a  source  (absorbing) 
region,  each  of  significant  optical  thickness.  While  most  methods  performed  equally  well 
in  the  source  region  with  an  incident  current,  fluxes  and  currents  became  less  accurate 
with  increasing  penetration  depth.  Of  all  of  the  methods,  the  EC  scheme  yielded  the  most 
accurate  solution  over  the  coarsest  meshes  (to  within  a  few  percent  of  the  reference 
solution).  In  addition  to  using  fewer  than  1/3  as  many  cells  as  the  linear  adaptive  method 
(the  closest  competitor),  the  EC  scheme  required  only  1/3  of  the  execution  time  required 
for  LA  to  obtain  an  equivalent  solution  using  a  coarse  mesh.  Again,  these  slab  geometry 
results  suggest  that  nearly  a  tenfold  savings  in  memory  storage  and  computational  cost  is 
possible  if  the  EC  scheme  is  implemented  in  two  dimensional  Cartesian  geometry. 

The  EC  method  encountered  problems  for  very  flat  flux  profiles  when  Sx  — >  0,  (the 

exponential  constant  P  is  close  to  zero),  although  this  was  easily  solved  by  using  the  LC 
method  in  such  instances.  Under  these  circumstances,  the  EC  method  is  asymptotically 
equivalent  to  the  LC  method. 
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In  the  strongly  diffusive  region  of  test  case  3,  EC  solutions  compared  well  with 
solutions  for  LC  and  LA.  Overall,  the  solutions  for  EC  were  more  accurate.  These 
limited  results  suggest  that  the  EC  method  behaves  correctly  in  the  thin  cell  diffusion 
limit.  Further,  although  yielding  a  solution  close  to  that  provided  by  LC  and  LA,  the  EC 
method  typically  required  more  than  twice  the  computational  effort  (mainly  due  to  root 
solving  requirements).  Improved  efficiency  in  root  solving  would  decrease  the 
disadvantage  EC  demonstrates  here,  while  increasing  the  advantage  of  using  EC  in  the 
other  test  cases. 

The  advantages  of  using  successive  scatters  instead  of  a  fixed  point  source  iteration 
(with  the  exponential  characteristic  scheme)  were  not  significant.  The  greatest  benefit 
occurred  as  optical  thicknesses  became  very  thin  in  absorbing  regions.  In  general, 
iteration  of  the  source  using  successive  scatters  performed  worse  than  fixed  point  over 
coarse  meshes  and  better  than  fixed  point  over  fine  meshes.  It  is  likely  that  the 
attenuation  of  successive  flignts  over  very  coarse  meshes  leads  to  increased  truncation 
error.  Since  accurate  calculations  are  desired  over  very  coarse  cells,  use  of  successive 
scatters  is  not  advantageous. 

Overall,  the  exponential  characteristic  method  proved  to  be  a  robust  scheme, 
performing  best  in  deep  penetration  problems,  when  |  SJ  »  SA.  Further,  considering  all 
spatial  quadrature  methods  used  in  this  research,  the  exponential  characteristic  scheme 
typically  yielded  the  most  accurate  solutions  over  the  coarsest  meshes  in  most  problems. 
This  was  especially  true  in  problems  having  several  regions  with  various  scattering  and 
absorption  properties.  In  most  cases,  the  exponential  characteristic  method  approached 
fourth  order  convergence  as  the  mesh  was  refined.  Thus,  numerical  diffusion  using  the 
EC  scheme  was  minor.  Even  though  the  EC  scheme  is  non-linear,  no  difficulties  in 
obtaining  convergence  were  encountered  for  any  problem  tested. 
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There  are  several  issues  concerning  exponential  characteristic  spatial  quadrature  that 
warrant  more  attention.  By  far,  further  optimization  of  the  root  solving  task  offers  the 
single  largest  improvement  in  performance  for  the  EC  method,  potentially  reducing  total 
execution  time  required.  Also,  better  treatment  of  root  solving  difficulties  when  P  values 
are  close  to  zero  is  warranted,  possibly  by  using  series  expansions  in  (3  in  such  cases. 
Moreover,  further  investigation  of  the  performance  of  the  EC  scheme  in  the  diffusion 
limit  is  warranted.  Finally,  the  utility  and  feasibility  of  conserving  second  spatial 
moments  for  this  scheme  might  be  investigated,  using  combinations  of  exponentials 
(hyperbolic  forms)  for  the  source  function.  In  any  event,  this  research  has  demonstrated 
that  the  exponential  characteristic  scheme  is  a  powerful  technique  for  solving  discrete 
ordinates  problems  in  slab  geometry.  Ultimately,  the  exponential  characteristic  scheme 
may  be  used  to  evaluate  transport  problems  once  thought  too  large  or  too  difficult  to  be 
adequately  solved  on  conventional  computer  systems. 
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All  slab  geometry  discrete  ordinates  transport  problems  were  solved  using  S8 

Gauss-Legendre  angular  quadrature.  Quadrature  direction  cosines  and  weights  are  given 
in  Table  6,  with  a  graphical  representation  in  Figure  20. 


Table  6 

Eight  Point  Single  Range  Gauss-Legendre  Quadrature  Constants 


m 

=  em° 

^ mRad 

\lm  =  Cos(6m) 

n 

163.8 

0.905  k 

-.960289856497536 

.101228536290376 

2 

142.8 

0.793  k 

-.796666477413627 

.222381034453374 

fl 

121.7 

0.676  n 

-.525532409916329 

.313706645877877 

D 

100.6 

0.559  n 

-.18343464249565 

.362683783378362 

5 

79.4 

0.441  Tt 

+.18343464249565 

.362683783378362 

6 

58.3 

0.324  7t 

+.525532409916329 

.313706645877877 

7 

37.2 

0.207  7t 

+.796666477413627 

.222381034453374 

8 

16.2 

0.090  7t 

+.960289856497536 

.101228536290376 
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Descriptions  of  three  types  of  boundaries  and  three  kinds  of  incident  currents 
available  in  the  multi-method  fixed  point  and  successive  scatters  discrete  ordinates  codes 
are  provided  here.  More  detailed  descriptions  of  each  are  available  in  the  literature 
(Lewis  and  Miller,  1984).  Boundary  conditions  made  available  in  the  discrete  ordinates 
codes  SN1D-ALL  and  SS1D-ALL  are  a  vacuum,  a  symmetric  albedo,  and  a  grey  albedo. 
Incident  left  or  right  boundary  currents  available  are  a  Lambertian  current,  an  isotropic 
surface  source  current,  or  a  collimated  beam  current. 

i.  Vacuum  Boundary 

If  a  vacuum  boundary  is  specified,  this  means  that  the  area  beyond  the  last  physical 
region  defined  for  the  problem  is  a  vacuum.  No  neutrons  can  be  scattered  or  reflected 
back  into  the  problem  after  crossing  the  physical  boundary.  This  is  often  called  a  "zero 
return  current"  condition. 

ii.  Symmetric  Albedo  Boundary 

At  a  symmetric  albedo  boundary,  for  each  amount  of  flux  leaving  across  this 
boundary,  an  albedo  a  (a  fraction  between  0  and  1)  of  the  original  outbound  flux  enters  in 
a  direction  corresponding  to  a  spectral  reflection.  A  symmetry  boundary  is  a  special  case 
of  the  symmetric  albedo  boundary,  where  the  albedo  a  =  1.0;  this  is  often  specified  when 
a  completely  symmetrical  problem  is  specified  and  a  solution  for  only  half  of  the  entire 
problem  is  necessary. 
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iii.  Grev  Albedo  Boundary 

The  grey  albedo  boundary  holds  that  an  albedo  a  times  the  outbound  flux  of  neutrons 
passing  across  the  surface  of  the  boundary  return  back  into  the  boundary  in  an  isotropic 
distribution,  uniformly  dispersed  over  all  angles.  This  is  accomplished  in  slab  geometry 
by  noting,  for  instance,  at  a  left  boundary,  the  entering  (positive)  current  is  given  by 

r(0)=l/op¥(°,4)f  (84) 

If  the  incoming  flux  is  a  constant  value  \)/+(0)  =  A  (as  specified  by  the  isotropic  condition) 
for  directions  0  <  |i  <  1,  then  the  positive  current  reduces  to 

(85) 

From  continuity  of  current  at  the  boundary,  which  yields 

y+(0)=/4^4(xT(0)  (86) 

This  states  that  for  a  grey  albedo  left  boundary,  the  outgoing  (negative)  current  is 
multiplied  by  four  times  the  albedo  factor  a  to  yield  an  isotropic  incoming  angular  flux 
over  0  <  p  <  1. 

A  special  case  of  the  grey  albedo  boundary  is  a  white  boundary,  where  a  =  1,  and  all 
exiting  current  returns  in  an  isotropic  distribution. 
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A  Lambertian  incident  current  is  due  to  particle  emissions  from  an  outside  isotropic 
source,  so  that  the  angular  flux  which  contributes  to  a  Lambertian  current  is  a  constant  F. 
Using  a  left  boundary  and  integrating  over  positive  angles,  the  following  equation  yields 
the  incident  current : 


•Ct(0)= 


F_ 

4 


Thus,  the  incident  current  attributable  to  each  direction  cosine  is 


(87) 


J,J  H)  =  I ^ncL 

A  similar  treatment  is  made  at  the  right  boundary. 


(88) 


v.  Isotropic  Surface  Source  Current 

An  isotropic  surface  source  is  a  planar  source  of  isotrop'c  emitters,  where  the  product 
of  the  angular  flux  and  direction  cosine  is  held  constant.  Therefore,  when  the  direction 
cosine  (i  is  close  to  zero,  the  angular  flux  is  highest  due  to  greatest  contribution  from  the 
isotropic  surface  sources.  Subsequently,  the  incident  current  attributable  to  each 
direction  cosine  at  a  left  boundary  is 

y,..(n)  =  ic,  m 


vi.  Collimated  Beam  Current 

Discrete  ordinates  methods  often  do  not  accurately  deal  with  highly  anisotropic  flux 


distributions,  e.g.  a  collimated  beam.  This  is  because  the  ordinates  of  a  standard 


quadrature  set,  symmetric  from  -1  to  1,  do  not  individually  provide  an  accurate  account 
of  the  forward  angular  flux  from  an  incident  current.  Under  these  circumstances,  the  best 
method  of  treating  a  collimated  beam  incident  in  a  direction  cosine  |i'  that  is  not  part  of 
the  quadrature  set  is  to  interpolate  the  incident  direction  between  applicable  directions  in 
the  quadrature  set.  Thus,  the  incident  current  in  ji'  is  shared  between  |i,  and  (l,  +  1,  where 
quadrature  weights  and  direction  cosines  are  used  to  adjust  the  magnitude  of  the  incident 
current  to  yield  accurate  quadratures. 
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Presented  below  are  the  test  cases  and  their  reference  solutions  using  fixed  point 
source  iteration  (solutions  obtained  using  successive  scatters  are  similar).  These 
reference  solutions  were  used  to  illustrate  the  performance  of  the  exponential 
characteristic  method  in  Section  VI.  Only  test  case  1  was  acquired  using  an  8  MHz  8086 
personal  computer  with  8087  math  support.  All  other  problems  were  solved  using  a  33 
MHz  80486  ISA  machine.  (The  80486  machine  tested  approximately  twenty  times  faster 
than  the  8086  machine). 


i.  Test  Case  1:  Thick  Absorber 


The  input  parameters  and  reference  solution  for  the  single  region  thick  absorber 
problem  are  given  below: 

Problem  file:  A.TMA1.IPT  Output  file:  B:TMA1-256.0UT 

Ident$:  MdlP  Mathews  90 
number  of  regions  =  1 
left  bdy  position  =  0 
type  ot  left  bdy  =  0  vacuum 
current  incident  at  left  boundary  =  1 
type  of  current  incident  at  left  bdy  =  0  isotropic  surface  Src 
region#  cR  SigmaR  SourceR  nc  Right  Bdy 
1  0.1000  1.0000D+00  O.OOOOD+OO  256.  16.0000 

type  of  right  bdy  =  0  vacuum 
current  incident  at  right  boundary  =  0 
type  of  current  incident  at  right  bdy  =-l  Lambertian 


Problem  file  :  A:TMA1.IPT  nk  Angular  Ordinates:  8 
Output  file  :  B:TMA  1-256. OUT  Quadrature  Method  :  S8 
Tk  file  :  B:TMA1-256.TKD  Negative  Flux  Fixups:  NO 
Max  Iterations  :  150  Reset  Old  Fluxes  :  YES 

File  Identifier:  MdlP  Mathews  90  Solution  Tolerance  :  1.00D-05 


00:45:20  1 1-01-1991  SN1D-ALL.BAS  Version  2.00  -  30  Oct  1991 

Linear  Characteristic  ->  LC  Source  Rotation  (Sx<=Sa)  Fixup  ALWAYS  Enabled 
Execution  Time  :  2.4488  min 

Converged  After  8  iterations,  MaxChangeObs  =  9.686617 14 1449785D-06 

J  plus  J  minus  J  net  Bdy  Flux 
Reg#  flux  ave  flux  x-moment 

+1.000000D+00  +2.633397D-02  +9.736660D-01  +3.022184D+00 
1  +6.761570D-02  -1.892749D-01 

+8.612 159D-09  +0.000000D+00  +8.6121 59D-09  +9.517948D-09 

Mom  Bals:  0th:  CONSERVED  MaxRelErr:  4.4702D-13<  1.0000D-10 
1st:  VIOLATION  1.4321D-10<  1.0000D-10 

Region  Bals:  CONSERVED  3.6304D-08<  1  .OOOOD-05 


ii.  Test  Case  2:  Moderate  Scatter  and  Source 


The  input  parameters  and  reference  solution  for  the  moderate  scatter  and 
absorber- source  problem  are  given  below: 


Problem  file:  D.TGSI.IPT  Output  file:  D:TGS  1 -256.0UT 


IdentS.  Ges  16mfp  Abs+Src 

number  of  regions  =  2 

left  bdy  position  =  0 

type  of  left  bdy  =  0  vacuum 

current  incident  at  left  boundary  =  0 

type  of  current  incident  at  left  bdy  =  0  isotropic  surface  Src 

region#  cR  SigmaR  SourceR  nc  Right  Bdy 

1  0.5000  1.0000D+00  0.0000D+O0  256.  16.0000 

2  0.1000  1 . 0000 D +00  5.0000D-01  256.  32.0000 
tyoe  of  right  bdy  =  0  vacuum 

current  incident  at  right  boundary  =  1 

type  of  current  incident  at  right  bdy  =  0  isotropic  surface  Src 
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IdentS:  Ges  16mfp  Abs+Diff 

number  of  regions  =  2 

left  bdy  position  =  0 

type  of  left  bdy  =  0  vacuum 

current  incident  at  left  boundary  =  0 

type  of  current  incident  at  left  bdy  =  0  isotropic  surface  Src 

reg  #  cR  SigmaR  SourceR  nc  Right  Bdy 

1  0  1000  1 . 0000 D +00  0.0000D+00  256.  16.0000 

2  0.9500  1.0000D+00  1.0000D+00  256.  32.0000 
type  of  right  bdy  =  0  vacuum 

current  incident  at  right  boundary  =  1 

type  of  current  incident  at  right  bdy  =  0  isotropic  surface  Src 
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Problem  file  :  D.TGS2.IPT  nk  Angular  Ordinates:  8 
Output  file  :  D:TGS2-256.0UT  Quadrature  Method  :  S8 
Tk  file  :  D:TGS2-256.TKD  Negative  Flux  Fixups:  NO 
Max  Iterations  :  200  Reset  Old  Fluxes  :  YES 

File  Identifier:  Gesl6mfp  Abs+Dif  Solution  Tolerance  :  1.00D-05 

03:55:40  12-02-1991  SN1D-ALL.BAS  Version  2.00  -  30  Oct  1991 

Linear  Characteristic  ->  LC  Source  Rotation  (Sx<=Sa)  Fixup  ALWAYS  Enabled 
Execution  Time  :  6.5388  min 

Converged  After  144  iterations,  MaxChangeObs  =  9.652863571 151075D-06 

J  plus  J  minus  J  net  Bdy  Flux 
Reg  #  flux  ave  flux  x-moment 

+0.000000D+O0  +3.550017D-08  -3.550017D-08  +3.913293D-08 

1  + 1 .38666 1 D-0 1  +3.786462D-0 1 

+4.43064 ID-02  +2.041 100D+00  - 1 .996793D+00  +3.803869D+00 

2  +1.544256D+01  +9.081632D-01 

+2.647280D+00  + 1 .000000D+00  +1.647280D+00  +7.970842D+00 

Mom  Bals:  0th:  CONSERVED  MaxRelErr:  8.9809D-13<  1.0000D-10 
1  st:  VIOLATION  4.4807D-07<  1 .0000D- 1 0 

Region  Bals:  CONSERVED  7.5017D-06<  1  .OOOOD-05 


Appendix  D;  Source  Code 


The  following  is  the  source  listing  for  the  multi-method,  slab  geometry,  fixed  point 
source  iteration  SN1D-ALL  code,  written  in  MicroSoft  QuickBASIC  4.5.  Except  for  the 
exponential  characteristic  scheme  implemented  by  the  author,  all  methods  were  derived 
from  individual  codes  originally  written  by  LCDR  K.  Mathews,  Ph.D.  at  the  Department 
of  Engineering  Physics,  Air  Force  Institute  of  Technology,  Wright  Patterson  AFB,  OH, 
45324.  As  discussed,  the  author  incorporated  all  codes  into  this  single  compilation  for  a 
uniform  analysis  approach.  The  successive  scatter  transport  code,  SS1D-ALL,  and  the 
error  compilation  code,  PROCSN02,  are  archived  with  LCDR  Mathews  at  the  above 
address. 


DECLARE  SUB  SetSwitchECtoLC  () 

DECLARE  SUB  StepDD  ( Jin* , Jout * , Sx* , dx* , Sigma* , mu* - Fa* ) 

DECLARE  SUB  StepSC  ( Sigma* , DeltaX* , mu* , FI* , Sa*  ,  Sx*  ,  Fr*  ,  Fa*  ,  Fx» ) 
DECLARE  SUB  StepLC  ( Sigma* , DeltaX# , mu* , FI* , Sa* , Sx*  ,  Fr*  ,  Fa*  ,  Fx» ) 
DECLARE  SUB  StepLN  (Sigma# , DeltaX# , mu# , FI* , Sa*  ,  Sx*  ,  Fr*  ,  Fa* , Fx* ) 
DECLARE  SUB  StepSA  ( Sigma* , DeltaX* , mu* , FI* , Sa* , Sx# , Fr* , Fa* , Fx* ) 
DECLARE  SUB  StepLA  ( S igma # , DeltaX# , mu* , FI# , Sa# , Sx* , Fr*  ,  Fa* , Fx* ) 
DECLARE  SUB  StepEC  ( Sigma* , DeltaX# , mu# , FI# , Sa# , Sx# , Fr* , Fa* , Fx* ) 
DECLARE  SUB  McmOStepCheck  ( FI* , Fr * , Fa* , eps# , Sa* , DeltaX* , mu# ) 
DECLARE  SUB  MomlStepCheck  ( Fl« , Fr * , Fa* , Fx* , eps* , Sx* , Del t aX* , mu* ) 
DECLARE  SUB  PrtTK  (labl$,  vrbl«(),  istart%,  n%) 

DECLARE  SUB  BeqinSets  () 

DECLARE  SUB  AltMenu  () 

DECLARE  SUB  makeP  (p*(),  e* ) 

DECLARE  FUNCTION  MIN*  (x«,  y«> 

DECLARE  FUNCTION  MAX*  (x“,  y« ) 

DECLARE  FUNCTION  Chocsel  (promptS) 

DECLARE  FUNCTION  ParseNthWordS  (a$,  n%) 

DECLARE  FUNCTION  Gx»  (x*,  ro«) 

DECLARE  FUNCTION  DdxG*  (x« ) 

CONST  False  =  0 
CONST  True  =  NOT  False 

DEFINT  I,  K,  N 

DEFDBL  A-G,  J,  L-M,  O-Z 

pi  =  4  *  ATN ( 1 ! ) 

'PROGRAM  SN ID- ALL. BAS 

'Discrete  Ordinates  (SN)  ID  Slab  Geometry  Program 

'for  research  at  the  Air  Force  Institute  of  Technology 

'by  Kirk  A.  Mathews,  Ph.D,  modified  by  Glenn  E.  Sjoden,  P.E. 

V$  =  "SN1D-ALL . BAS  Version  2.10  -  10  Dec  1991“ 
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Spatial  Quadrature  Methods  Supported  in  this  Code  -- 
'  DD-Diamond  Difference  (from  V-2.05,  19  Feb  90) 

'  LD-Linear  Discontinuous  (from  V-2.04,  4  Feb  90) 

'  SC-Step  Characterist ic  (from  V-2.04,  1  Mar  90) 

'  LC-Linear  Character ist ic  (Modified  from  LA  below) 

'  LN-Linear  Nodal  (from  V-2.04,  4  Feb  90) 

'  SA-Step  Adaptive  (from  V-2.04,  4  Feb  90) 

'  LA-Linear  Adaptive  (from  V-2.04,  4  Feb  90) 

'  -  by  Kirk  A.  Mathews,  Ph.D. 

'  EC-Exponential  Characteristic  (V-1.20,  16  Oct  91) 

'  -  by  Glenn  E.  Sjoden,  P.E. 

WIDTH  LFRINT  80 

'Arrays  used  in  defining  the  problem  -- 
'  Xbdy(ir)  is  right  bdy  of  region  ir 

'  cR(ir)  is  scatter  to  total  cross-section  ratio  in  region  ir 
'  SigmaR(ir)  is  total  cross-sect  ion  in  region  ir 
'  source(ir)  is  volumetric  source  in  region  ir 
'  nc(ir)  is  number  of  cells  in  region  ir 

'  mu(k)  is  k'th  "discrete  ordinate"  (direction  in  quadrature  set) 

'  w(k)  is  weight  associated  with  mu(k)  in  angular  quadrature  set 
'  dx(ix)  is  thickness  of  cell  ix 

'  Fa(ix,k)  is  flux  (psi)  in  directio,.  mu(k)  averaged  over  cell  ix 

'  Fx(ix,k)  is  first  moment  of  psi  in  direction  mu(k)  over  cell  ix 

'  J(ix,k)  is  current  (mu*psi)  in  direction  mu(k) 

'  through  right  side  of  cell  ix 

'  Note:  j  is  posit ive  for  positive  mu  and  neg .  for  neg .  mu 

'  FluxA(ix)  is  scalar  flux  averaged  over  cell  ix 
'  FluxX(ix)  is  first  moment  of  scalar  flux  over  cell  ix 

'  sigma(ix)  is  total  cross-section  in  cell  ix 

c(ix)  is  scatter  to  total  cross-section  ratio  in  cell  ix 
'  source (xx)  is  volumetric  source  in  cell  ix 

'  Sa(ix)  is  average  source  term  (s  =  c*f lux+source/sigma)  in  cell  ix 
'  Sx(ix)  is  change  in  Sa  across  cell  ix 
'  lprint  ix,  nx 

Problemfile$  =  Par seNthWord (( COMMANDS )  ,  1) 

'returns  first  parameter  found  on  command  line 
DO 

GOSUB  CommandProf i le 
SELECT  CASE  CmdProfS 
CASE  "S" 

GOSUB  OpenProblemFi le 
NewProblem%  =  True 
SwECt  cLCSet  Pt  =  .000012 
SwitchEC%  =  False 
GOSUB  ReadProblemData 
CLOSE  “1  'close  input  file 
DevEcho$  =  "SCREEN” 

DO 

CLS 

GOSUE  EchoProblem 
GOSUB  Adiust ProblemData 
LOOP  UNTIL  EnterNC%  =  False 

ff%  =  Choose ("Use  Form  Feeds  in  Summaries?  (Y/N) :  ") 

DevEcho$  =  "LPT1" 

GOSUB  Va 1 idat eProblem 
PRINT 

IF  Choose ( " Echo  Input  Data?  (Y/N):  ")  THEN 
DevEchoS  =  " LPT1 " 

GOSUB  EchoProblem 
END  IF 
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GOSUB  ResetMomerrs 
DO 

GOSUB  InputMethodParameters 
GOSUB  SetClocks 
DevEcho$  =  "LPT1" 

GOSUB  EchoMethod 
DevEcho$  =  "SCREEN" 

GOSUB  EchoMethod 

GOSUB  InitializeProblem 

NewProblem%  =  False 

GOSUB  SnSolve 

GOSUB  CalculateResults 

GOSUB  Ver i fyRegionsByBalanceEquat ion 

StopC lock  =  TIMER 

ExecMin  =  (StopClock  -  StartClock)  /  60! 

DevEcho$  =  “LPT1" 

GOSUB  LPr intRegionSummary 
GOSUB  CreateTKDataFi le 
GOSUB  ResetMomerrs 
SwitchEC%  =  False 

LOOP  WHILE  Choose ( "Solve  sam;  problem  with  different  Sn?  (y/n)  : 

"  ) 

PRINT 

ProblemfileS  =  "" 

CASE  “A” 

ff%  =  False 

CALL  BeginSets 

CALL  AltMenu 

NewProblem%  =  True 

GOSUB  ReadProblemData 

CLOSE  «  1  'close  input  -lie 

DevEchoS  =  "SCREEN" 

DO 

CLS 

GOSUB  EchoProblem 
GOSUB  AdjustProblemData 
LOOP  UNTIL  EnterNCt  =  False 
DevEchoS  =  "SCREEN" 

GOSUB  ValidateProblem 
IF  UCASES (Qoutf $)  =  "Y”  THEN 
DevEchoS  =  "OUT" 

GOSUB  EchoProblem 
END  IF 

GOSUB  ResetAgain 
Re  f luxA : 

DO 

IF  aga in%  THEN 
CALL  AltMenu 
DevEchoS  =  "SCREEN" 

DO 

CLS 

GOSUB  EchoProblem 
GOSUB  AdjustProblemData 
LOOP  UNTIL  EnterNC%  =  False 
DevEchoS  =  "SCREEN" 

GOSUF  ValidateProblem 
IF  UCASES (Qoutf $)  =  "Y"  THEN 
DevEchoS  =  "OUT" 

GOSUB  EchoProblem 
END  IF 

GOSUB  ResetAgain 
END  IF 


86 


GOSUB  Ir.putMethodParameters 
GOSUB  SetClocks 

DevEcho$  =  " SCREEN" :  Eout $  =  “ FULL" 

GOSUB  EchoMethod 

LOCATE  12,  12:  PRINT  “Processing  ->  TotDifMeth$;  "  <-“ 

IF  UCASE$ (Qoutf $)  =  “Y"  THEN 
DevEcho$  =  " OUT" 

IF  itim  =  0  THEN 
Eout  $  =  “FULL" 

ELSE 

Eout  $  =  “ PARTIAL" 

END  IF 

GOSUB  Echowethod 
END  IF 

GOSUB  InitializeProblem 
IF  ResetAl 1 Fluxes $  =  “YES"  THEN 
NewProblem%  =  True 
ELSE 

NewProblem%  =  False 
END  IF 

GOSUB  SnSolve 

GOSUB  CalculateResults 

GOSUB  Ver i f yRegionsByEalanceEquot ion 

StopClock  =  TIMER 

ExecMin  =  (StopClock  -  StartClock)  /  60! 

GOSUB  Incremitim 

DevEchoS  =  "SCREEN":  Eout$  =  “FULL" 

GOSUB  EchoMethod 
GOSUE  LPr intReaionSummary 
IF  UCASES (QoutfS)  =  “ Y "  THEN 
DevEchoS  =  "OUT" 

GOSUB  LPr int Reg ionSummary 
END  IF 

GOSUB  Creat eTKDataFi le 
GOSUB  ToqgleDlfMeth 
GOSUB  ResetMomerrs 

LOOP  WHILE  INSTR ( " DDLDSCLCLNSALAEC " ,  DlfMethS)  AND  LEN ( Di f Meth$ )  =  2 
PpvErho?  =  "SCREEN” 

GOSUB  RunSummary 
IF  UCASES ( Qout  f  $ )  =  “Y"  THEN 
DevEchoS  =  "OUT" 

GOSUB  RunSummary 
END  IF 
BEEP:  BEEP 

IF  Choose i "Solve  same  problem  with  different  Sn?  (y/n):  ” )  THEN 
aaain%  =  True 
GOTO  RefluxA 
END  IF 
END  SELECT 

LOCF  WHILE  Choose i "Solve  Another  Problem?  (y/n):  ") 

END 


CommandFrof i le : 

CLS  :  PRINT  VS:  PRINT 
DO 

INPUT  "Type  of  Command  Profile  (S-Standard/  A-Alternate)  ",  CmdProfS 
CmdProf $  =  UCASES ( LEFTS ( LTRIM$ (CmdProfS )  ,  1)) 

LOOP  UNTIL  INSTR ( “ SA “ ,  CmdProfS)  AND  LEN ( CmdProf $ )  =  1 
RETURN 
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SetClocks : 

StartDate$  =  DATE$ :  StartTime$  =  TIME$ 

StartClock  =  TIMER:  CLS 
RETURN 

ResetAgain : 
again%  =  False 
SwrtchEC%  =  False 

REDIM  Rtim ( 1  TO  8),  RelRtim(l  TO  8) 

REDIM  RErranal ( 1  TO  8) ,  RelRErranal ( 1  TO  8) 
itim  =  0 
RETURN 

Incremitim: 
itim  =  itim  +  1 
Rtim (itim)  =  ExecMin 
RErranal ( itim)  =  MaxRegrerr 
CLS 

RETURN 

ResetMomerrs : 

MomOIntFlag$  =  -  Oth:  CONSERVED":  MaxMomOrerr  =  0 
MomlIntFlag$  =  “  1st:  CONSERVED":  MaxMomlrerr  =  0 
RETURN 

ToggleDifMeth : 

ilotogl  =  ilotogl  +  3:  ihitogl  =  ihitogl  +  3 
DifMeth$  =  LTRIM$ (MID$ (TotDi fMeth$ ,  ilotogl,  ihitogl)) 
DifMeth$  =  UCASE$ (LEFT$ (LTRIM$ (DifMeth$) ,  2)) 

RETURN 

OpenProblemFile : 

CLS 

PRINT  V$ 

PRINT 

IF  Problemf ile$  =  THEN  INPUT  "File  to  read  for  input 
Problemf ile$ 

OPEN  Problemf ile$  FOR  INPUT  AS  #1 
'check  for  correct  file: 

INPUT  #1,  Ident$ 

PRINT  Ident$ 

DO 

INPUT  "Correct  File?  (Yes/No/Quit):  ",  a$ 
a$  =  UCASE$(LEFT$(LTRIM$(a$) ,  1)) 

LOOP  UNTIL  INSTR ( " YNQ" ,  a$)  AND  LEN (a$ )  =  1 
SELECT  CASE  a$ 

CASE  "N" 

CLOSE  #1 

Problemf ile$  =  "" 

GOTO  OpenProblemFile 
CASE  "Q" 

END 

CASE  "Y" 

CLS 

CASE  ELSE 
BEEP 

PRINT  "ERROR:  Unsupported  choice  in  OpenProblemFil 
STOP 
END  SELECT 
RETURN 


ReadProblemData : 

INPUT  #1,  nr  'number  of  regions 

REDIM  Xbdy (nr) ,  cR(nr),  SigmaR(nr),  SourceR(nr),  nc(nr) 

INPUT  #1,  Xbdy ( 0 )  'left  bdy  position 

INPUT  #1,  tlb  'type  of  left  bdy  — 

' 0  =  vacuum 

'(0,+l]  =  symmetric  albedo 
'{specular  reflection  factor) 

' [-1,0)  =  grey  albedo 
' {Lambertian  reflection  factor) 

INPUT  #1,  jinclb  'current  incident  at  left  boundary 

INPUT  #1,  tinclb  'type  of  current  inc  at  left  bdy  -- 

' -1  =  Lambertian 
'0  =  isotropic  surface  source 

'(0,+l)  =  abs(mu)  of  collimated  incident  beam 
FOR  ir  =  1  TO  nr 

INPUT  #1,  cR(ir),  SigmaR ( ir ) ,  SourceR(ir),  nc(ir),  Xbdy(ir) 

NEXT  ir 

INPUT  #1,  trb,  jincrb,  tincrb 
RETURN 

AdjustProblemData : 

EnterNC%  =  Choose ( "Manually  enter  values  for  #  cells  /  region?  (y/n) 
") 

IF  EnterNC%  THEN 
FOR  ir  =  1  TO  nr 

PRINT  ”  ir;  TAB (10); 

PRINT  USING  “##.####";  cR(ir); 

PRINT  TAB (20) ; 

PRINT  USING  “##.####AAAA";  SigmaR(ir); 

PRINT  TAB (34) ; 

PRINT  USING  "##  .####AAAA“ ;  SourceR(ir); 

PRINT  TAB (48) ; 

PRINT  USING  ”->  ###.  nc(ir); 

PRINT  TAB ( 58 ) ; 

PRINT  USING  “###.####";  Xbdy(ir) 

IF  EnterNC%  THEN 

PRINT  “Region”;  ir; 

INPUT  "Number  of  cells?  nc  =  ",  nc(ir) 

PRINT 
END  IF 
NEXT  ir 
END  IF 
PRINT 
RETURN 

Validate Problem: 
bad%  =  False 
SELECT  CASE  DevEcho$ 

CASE  "LPT1 " 

OPEN  " LPT1 : "  FOR  OUTPUT  AS  #4 
CASE  "SCREEN" 

OPEN  " SCRN : “  FOR  OUTPUT  AS  #4 
CASE  "OUT" 

OPEN  Out f i le$  FOR  APPEND  AS  #4 
END  SELECT 
FOR  ir  =  1  TO  nr 

IF  Xbdy ( i r )  <=  Xbdy(ir  -  1)  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  " Xbdy ( " ;  ir;  ")  =  ";  Xbdy(ir) 

bad%  =  True 
END  IF 
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IF  cR(ir)  <  0  OR  cR(ir)  >  1  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  " cR ( " ;  ir;  ")  =  * ;  cR(ir) 
bad%  =  True 
END  IF 

IF  SigmaR(ir)  <=  0  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  SigmaRC;  ir;  ")  =";  SigmaR(ir) 
bad%  =  True 
END  IF 

IF  SourceR ( ir )  <  0  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  SourceRC;  ir;  ")  =";  SourceR  (ir) 
bad%  =  True 
END  IF 

IF  nc(ir)  <=  0  THEN 

PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  nc ( " ;  ir;  ")  =*;  nc{ir) 
bad%  =  True 
END  IF 
NEXT  ir 

IF  ABS(tlb)  >  1  THEN 

PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  type  of  left  boundary  --  tlb  =";  tlb 
bad%  =  True 
END  IF 

IF  ABS(trb)  >  1  THEN 

PRINT  #4,  "Bad  Input:" 

PRINT  #4,  “  type  of  right  boundary  --  trb  =";  trb 

bad%  =  True 
END  IF 

IF  jinclb  <  0  THEN 

PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  negative  incident  current  at  left  --  jinclb  =" ;  jinclb 
bad%  =  True 
END  IF 

IF  jincrb  <  0  THEN 

PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  negative  incident  current  at  right  --  jincrb  =";  jincrb 
bad%  =  True 
END  IF 

IF  tinclb  <  -1  OR  (tinclb  >  -1  AND  tinclb  <  0)  OR  tinclb  >  1  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  type  of  incident  current  at  left  --  tinclb  =";  tinclb 

bad%  =  True 
END  IF 

IF  tincrb  <  -1  OR  (tincrb  >  -1  AND  tincrb  <  0)  OR  tincrb  >  1  THEN 
PRINT  #4,  "Bad  Input:" 

PRINT  #4,  "  type  of  incident  current  at  right  --  tincrb  =";  tincrb 

bad%  =  True 
END  IF 

IF  bad%  =  True  AND  (DevEcho$  =  “LPT1")  THEN 
PRINT  #4,  CHR$ ( 12 ) 

CLOSE  #4 
END 

ELSEIF  bad%  =  True  THEN 
CLOSE  #4 
END 
END  IF 
CLOSE  *4 
RETURN 
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EchoProblem: 

SELECT  CASE  DevEcho$ 

CASE  "LPT1" 

OPEN  " LPT1 : “  FOR  OUTPUT  AS  #4 
CASE  "SCREEN" 

OPEN  " SCRN : "  FOR  OUTPUT  AS  #4 
CASE  "OUT" 

OPEN  Outf ile$  FOR  APPEND  AS  #4 
END  SELECT 

PRINT  #4,  ******************************************************** 

*****★***★★**★» 

PRINT  #4,  "Problem  file:  UCASE$ ( Problemf ile$ ) ;  TAB(44);  "Output 

file:  UCASE$ (Outf ile$ ) 

PRINT  #4,  ******************************************************** 
PRINT  #4, 

PRINT  #4,  " Ident$ :  ";  Ident$ 

PRINT  #4,  "number  of  regions  =";  nr 
PRINT  #4,  "left  bdy  position  =";  Xbdy(O) 

PRINT  #4,  "type  of  left  bdy  =  ";  tlb; 

SELECT  CASE  tlb 
CASE  0 

PRINT  #4,  "  vacuum" 

CASE  0  TO  1 

PRINT  #4,  "  symmetric  albedo  {specular  reflection  factor}" 

CASE  -1  TO  0 

PRINT  #4,  "  grey  albedo  {Lambertian  reflection  factor}" 

CASE  ELSE 

PRINT  #4,  "  invalid" 

END  SELECT 

PRINT  #4,  "current  incident  at  left  boundary  =  ";  jinclb 
PRINT  #4,  "type  of  current  incident  at  left  bdy  =";  tinclb; 

SELECT  CASE  tinclb 
CASE  -1 

PRINT  #4,  "  Lambertian" 

CASE  0 

PRINT  #4,  “  isotropic  surface  Src” 

CASE  0  TO  1 

PRINT  #4,  B  abs(rmu)  of  collimated  incident  beam” 

CASE  ELSE 

PRINT  #4,  ■  invalid" 

END  SELECT 

PRINT  #4,  "region  #";  TAB(13);  "cR";  TAB(23);  "SigmaR";  TAB(37); 
PRINT  #4,  "SourceR";  TAB(49);  *nc  ";  TAB(58);  "Right  Bdy" 

FOR  ir  =  1  TO  nr 

PRINT  #4,  "  ";  ir;  TAB(IO); 

PRINT  #4,  USING  ”##.####";  cR(ir); 

PRINT  *4,  TAB (20); 

PRINT  #4,  USING  "  ##  .  ####~~~/'  "  ;  SigmaR(ir); 

PRINT  #4,  TAB (34) ; 

PRINT  #4,  USING  "  ##  .  "  ;  SourceR(ir); 

PRINT  #4,  TAB (48) ; 

PRINT  #4,  USING  "###.";  nc(ir); 

PRINT  #4,  TAB ( 58 ) ; 

PRINT  #4,  USING  "###.####";  Xbdy(ir) 

NEXT  ir 

PRINT  #4,  "type  of  right  bdy  =";  trb; 

SELECT  CASE  TypR 
CASE  0 

PRINT  #4,  "  vacuum" 

CASE  0  TO  1 

PRINT  #4,  "  symmetric  albedo  {specular  reflection  factor}" 


CASE  -1  TO  0 

PRINT  #4,  “  grey  albedo  {Lambertian  reflection  factor}" 

CASE  ELSE 

PRINT  #4,  "  invalid" 

END  SELECT 

PRINT  #4,  "current  incident  at  right  boundary  =";  jincrb 
PRINT  #4,  "type  of  current  incident  at  right  bdy  =";  tincrb; 

SELECT  CASE  tincrb 
CASE  -1 

PRINT  #4,  "  Lambertian" 

CASE  0 

PRINT  #4,  *  isotropic  surface  Src" 

CASE  0  TO  1 

PRINT  #4,  "  abs(rmu)  of  collimated  incident  beam" 

CASE  ELSE 

PRINT  #4,  *  invalid* 

END  SELECT 

IF  ff%  THEN  PRINT  #4,  CHR$(12)  ELSE  PRINT  #4,  :  PRINT  #4, 

SELECT  CASE  DevEcho$ 

CASE  ” LPT1 " ,  "SCREEN",  "OUT" 

CLOSE  #4 
END  SELECT 
RETURN 

InputMethodParameters : 

SELECT  CASE  CmdProf$ 

CASE  “S" 

PRINT 

DO 

PRINT  "Select  Spatial  Quadrature  Method  :" 

PRINT 

PRINT  "  DD  -  Diamond  Difference* 

PRINT  “  LD  -  Linear  Discontinuous" 

PRINT  "  SC  -  Step  Characteristic" 

PRINT  "  LC  -  Linear  Characteristic" 

PRINT  "  LN  -  Linear  Nodal" 

PRINT  "  SA  -  Step  Adaptive" 

PRINT  “  LA  -  Linear  Adaptive" 

PRINT  "  EC  -  Exponential  Characteristic" 

PRINT 

INPUT  "Choree?  (DD, LD, SC, LC, LN, SA, LA, EC) :  ",  DifMeth$ 

DifMeth$  =  UCASE$ (LEFT$ (LTRIM$ (DifMeth$)  ,  2)) 

LOOP  UNTIL  INSTR ( " DDLDSCLCLNSALAEC "  ,  DifMeth$)  AND  LEN ( Di f Meth$ )  =  2 
SELECT  CASE  DifMeth$ 

CASE  "DD" 

PRINT  "For  Diamond  Difference  Spatial  Quadrature:" 
fixup%  =  Choose("Use  Negative  Flux  Fixup?  (y/n)  :  ") 

PRINT 
CASE  "LN" 

PRINT  "For  Linear  Nodal  Spatial  Quadrature:" 

fixup%  =  Choose("Use  Scalar  Flux  Rotation  Fixup?  (y/n):  ") 

PRINT 

CASE  "LD",  "SC",  "LC",  "SA",  “LA" 

CASE  "EC" 

PRINT 

CALL  SetSwitchECtoLC 
END  SELECT 
DO 

INPUT  "Max  number  of  iterations:  IterMax  =  ",  Itermax% 

LOOP  UNTIL  Itermax  >  0 
PRINT 
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PRINT  “Criterion  for  convergence  is  that  over  an  iteration, " 

PRINT  “the  relative  change  in  each  cell  scalar  flux  <=  Change" 

DO 

INPUT  “Convergence  Criterion:  (0  <  Change  <  1)  Change  =  ", 

Change 

LOOP  UNTIL  0  <  Change  AND  Change  <  1 

PRINT 

DO 

INPUT  “Tolerance  for  Moment  Balance  Comparisons  (0<  MomTol  <  1)“ 

MomTo 1 

LOOP  UNTIL  Change  >  0# 

PRINT 

DO 

PRINT  “Select  Angular  Quadrature  Method:* 

PRINT  “  S  -  Single  Range  Gauss-Legendre“ 

PRINT  “  D  -  Double  Range  Gauss-Legendre" 

PRINT  “  M  -  Composite  Midpoint  Rule" 

PRINT 

INPUT  “Choice?  (S/D/M):  “,  QuadMeth$ 

QuadMethS  =  UCASE$ (LEFTS (LTRIMS (QuadMethS ) ,  D) 

LOOP  UNTIL  INSTR ( " SDM“ ,  QuadMethS)  AND  LEN(QuadMethS)  =  1 
PRINT 
CASE  "A" 

END  SELECT 

SELECT  CASE  QuadMethS 
CASE  “S“ 

GOSUB  SingleRangeGauss 
CASE  "D" 

GOSUB  DoubleRangeGauss 
CASE  "M" 

GOSUB  Compos iteMidpointQuadrature 
CASE  ELSE 
BEEP 

PRINT  "ERROR:  Illegal  choice  in  InputMethodParams " 

STOP 
END  SELECT 
RETURN 

SingleRangeGauss : 

SELECT  CASE  CmdProfS 
CASE  "S" 

DO 

PRINT  "For  Single-Range  Gauss  Sn,  n  is  total  #  of  mu’s." 

PRINT  "Supported  Orders  are  n  =  2,  4,  6,  8,  10  and  12." 

PRINT 

INPUT  "Enter:  n  =  ",  nk 

LOOP  UNTIL  (nk  >  =  2)  AND  (nk  <=  12)  AND  (nk  MOD  2=0) 

PRINT 
CASE  "A" 

END  SELECT 
REDIM  w(nk)  ,  mu(nk) 
nkO  =  nk 
GOSUB  EvenGauss 
RETURN 

EvenGauss : 

SELECT  CASE  nkO 
CASE  2 

mud)  =  -.577350269189626#:  w(l)  =  1# 

CASE  4 

mud)  =  -.861136311594053#:  w(l)  =  .347854845137454# 
mu  ( 2 )  =  -.339981043584856#:  w(2)  =  .652145154862546# 
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CASE  6 
mu  ( 1 ) 
mu  ( 2 ) 
mu  ( 3 ) 
CASE  8 
mu  ( 1 ) 
mu  (2) 
mu  ( 3 ) 
mu  (4) 
CASE  10 
mu  ( 1 ) 
mu  (2 ) 
mu  ( 3 ) 
mu  (4 ) 
mu  ( 5 ) 
CASE  12 
mu  ( 1 ) 
mu  (2) 
mu  ( 3 ) 
mu  ( 4 ) 
mu  ( 5 ) 
mu  ( 6 ) 


.9324695142031521#:  w(l)  = 
.661209386466265#:  w<2)  = 
.238619186083197#:  w(3)  = 

.960289856497536#:  w{l)  = 
.7966664774136269#:  w(2)  = 
.525532409916329#:  w<3)  = 

.18343464249565#:  w(4)  = 

.973906528517172#:  w(l)  = 
.865063366688985#:  w(2)  = 
.679409568299024#:  w(3)  = 
.433395394129247#:  w(4)  = 
.148874338981631#:  w(5)  = 

.981560634246719#:  w(l)  = 
.904117256370475#:  w(2)  = 
.769902674194305#:  w(3)  = 
.587317954286617#:  w(4)  = 
.36783149899818#:  w(5)  = 
.125233408511469#:  w(6)  = 


.17132449237917# 

.360761573048139# 

.467913934572691# 

=  .101228536290376# 
=  .222381034453374# 
=  .313706645877877# 
=  .362683783378362# 

.066671344308688# 

.149451349150581# 

.219086362515982# 

.269266719309996# 

.295524224714753# 

.047175336386512# 

.106939325995318# 

.160078328543346# 

.203167426723066# 

.233492536538355# 

.249147045813403# 


CASE  ELSE 
BEEF 

PRINT  “ERROR:  Gauss -Legendre  with  NK  =  ";  nkO;  "  not  supported. 
STOP 
END  SELECT 

FOR  k  =  nkO  TO  nkO  \  2  +  1  STEP  -1 
mu(k)  =  -mu(nk0  -  k  +  1) 
w(k)  =  w(nk0  -  k  +  1) 

NEXT  k 
RETURN 


DoubleRangeGauss  : 

SELECT  CASE  CmdProf$ 

CASE  "S" 

DO 

PRINT  "For  Double-Range  Gauss  Sr.,  n  is  #  of  mu's  in  each  range 
PRINT  "Supported  Orders  are  n  =  1,  7,  3,  4,  6,  8,  10  and  12." 
INPUT  "Enter:  N  =  ",  nk 

LOOP  UNTIL  nk  >=  1  AND  nk  <=  12  AND  ( (nk  MOD  2=0)  OR  nk  <  4) 
CASE  "A" 

END  SELECT 
nkO  =  nk 
nk  =  2  *  nkO 
PRINT 

REDIM  w(nk),  mu(nk) 

SELECT  CASE  nkO 
CASE  1 

mu(l)  =  0:  w(l)  =  2 
CASE  3 

mu ( 1 )  =  -.7745966692#:  w(l)  =  .5555555556# 
mu ( 2 )  =  0:  w(2)  =  .8888888888000001# 
mu ( 3 )  =  -mu (1):  w(3)  =  w(l) 

CASE  2,  4,  6,  8,  10,  12 
GOSUB  EvenGauss 
CASE  ELSE 
BEEP 

PRINT  "ERROR:  Illegal  nkO  in  DoubleRangeGauss." 

STOP 
END  SELECT 

'Shift  and  scale  the  above  single-range  quad  set  into  the  interval 


1) :  w(k) 


-5  *  w(k) 


(-1,0) 

FOR  k  =  1  TO  nkO 

mu(k)  =  .5  *  (mu(k)  - 
NEXT  k 

'mirror  to  the  other  single-range  quad  set,  for  interval  (0,+l) 
FOR  k  =  nkO  +  1  TO  nk 

mu(k)  =  -mu(nk  +  1  -  k):  w(k)  =  wink  +  1  -  k) 

NEXT  k 
RETURN 

Compos iteMidpointQuadrature : 

SELECT  CASE  CmdProf$ 

CASE  “S" 

DO 

PRINT  “For  Composite  Midpoint  Rule  Sn,  n  is  total  #  of  mu's. 
PRINT  "The  mu's  are  at  the  center  of  equally  sized  intervals 
PRINT  "  and  are  given  equal  weights." 

PRINT  "Supported  Orders  are  n=2,  4,  6,  ..." 

INPUT  “Enter:  N  =  ",  nk 
LOOP  UNTIL  nk  >  0  AND  ( (nk  MOD  2)  =0) 

PRINT 
CASE  "A" 

END  SELECT 

REDIM  w(nk),  mu(nk) 

FOR  k  =  1  TO  nk 

mu(k)  =  -1  +  (2  *  k  -  1)  /nk 
w(k)  =  1  /  nk 
NEXT  k 
RETURN 

EchoMethod : 

SELECT  CASE  DevEcho$ 

CASE  " LPT1 “ 

OPEN  ” LPT1 : “  FOR  OUTPUT  AS  #4 
CASE  "SCREEN" 

OPEN  “ SCRN : “  FOR  OUTPUT  AS  #4 
CASE  "OUT" 

OPEN  Out f i le$  FOR  APPEND  AS  #4 
END  SELECT 

SELECT  CASE  DifMeth$ 

CASE  "DD” 

a2$  =  "Diamond  Difference" 

IF  f ixup%  THEN 

a3$  =  "  with  Negative  Flux  Fixup." 

ELSE 

a3$  =  "  with  NO  Flux  Fixup." 

END  IF 
CASE  "LD" 

a2$  =  "Linear  Discontinuous" 
a3$  =  ” " 

CASE  "SC" 

a2$  =  "Step  Characteristic" 
a3  $  =  "" 

CASE  " LC “ 

a2$  =  "Linear  Characteristic" 
a3  $  =  *" 

CASE  "LN" 

a2$  =  "Linear  Nodal" 

IF  f ixup%  THEN 

a3$  =  "  with  Scalar  Flux  Rotation  Fixup." 

ELSE 

a3$  =  "  with  NO  Flux  Rotation  Fixup." 
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END  IF 
CASE  "SA" 

a2$  =  "Step  Adaptive" 
a3  $  =  “ “ 

CASE  “LA" 

a2$  =  “Linear  Adaptive" 
a3  $  =  “" 

CASE  "ECH 

a2$  =  "Exponential  Characteristic" 
a3  $  =  *" 

END  SELECT 

SELECT  CASE  CmdProf$ 

CASE  "S" 

PRINT  #4,  ********************************************************** 

PRINT  *4,  StartTime$;  •  StartDate$;  "  " ;  V$ 

PRINT  #4,  "Solution  by  Discrete  Ordinates  (Sn)  with  n  =  ";  nk; 

PRINT  #4,  "mu's,  total," 

SELECT  CASE  QuadMeth$ 

CASE  "S" 

al$  =  "  with  Single-Range  Gauss-Legendre  Angular  Quadrature" 
CASE  "D" 

al$  =  "  with  Double-Range  Gauss-Legendre  Angular  Quadrature" 
CASE  "M" 

al$  =  "  with  Composite  Midpoint  Angular  Quadrature" 

CASE  ELSE 
BEEP 

PRINT  "ERROR:  Unsupported  angular  quadrature  (QuadMeth$)  in 
EchoMethod . “ 

STOP 
END  SELECT 
PRINT  #4,  al$ 

PRINT  #4,  LTRIM$(*and  "  +  a2$  +  “  Spatial  Quadtrature"  +  a3$  + 

PRINT  #4,  "Convergence  criterion  on  scalar  fluxes:  RelChange  < 

Change  =";  Change 
PRINT  *4, 

CASE  "A" 

IF  Eout $  =  "FULL"  THEN 

PRINT  #4,  "********************************************************* 
•» 

PRINT  #4,  "Problem  file  :  ";  UCASE$ ( Problemf ile$ ) ;  TAB(44);  "nk 
Angular  Ordinates  : " ;  nk 

PRINT  #4,  "Output  file  :  - ;  UCASE$ (Out f ile$ ) ;  TAB (44); 

"Quadrature  Method  :  ";  UCASE$ (QuadMeth$ )  +  LTRIM$ (STR$ (nk) ) 

PRINT  #4,  "Tk  file  :  ";  UCASE$ (TKf ile$ ) ;  TAB(44);  "Negative 

Flux  Fixups  :  ";  UCASE$ ( FI f ix$ ) 

PRINT  <*4,  “Max  Iterations  :";  Itermax%;  TAB  (44);  "Reset  Old  Fluxes 
:  ";  Reset A1 lFluxesS 

PRINT  *4,  "File  Identifier  :  ";  Ident$;  TAB(44);  "Solution  Tolerance 
PRINT  *4 ,  USING  Change 

PRINT  #4,  ********************************************************** 

PRINT  #4,  StartTime$ ;  TAB(17);  StartDate$;  TAB(35);  V$ 

PRINT  #4, 

IF  DevEcho$  =  "OUT"  THEN  PRINT  #4,  LTRIM$(a2$  +  "  ->  "); 

ELSE 

PRINT  <*4,  >********************************************************* 

****************** 

PRINT  «4,  StartTimeS;  TAB (17);  StartDate$;  TAB (35);  V$ 

PRINT  *»4  , 

PRINT  *4,  LTRIM$(a2$  +  “  ->  "); 


96 


END  IF 
END  SELECT 
CLOSE  #4 
RETURN 

InitializeProblem: 

IF  NewProblem%  THEN 
nx  =  0 

FOR  ir  =  1  TO  nr 
nx  =  nx  +  nc(ir) 

NEXT  ir 

REDIM  Fa(nx,  nk) ,  Fx(nx,  nk) ,  J(nx,  nk) ,  Jinc(nk) 

REDIM  FluxA(nx),  FluxX(nx),  Sigma(nx),  c(nx),  Source(nx) 

REDIM  Sa(nx),  Sx(nx),  dx(nx),  CellCtrX(nx) 
ix  =  0 

FOR  ir  =  1  TO  nr 

WidthX  =  (Xbdy(ir)  -  Xbdy(ir  -  1))  /  nc(ir) 

FOR  ic  =  1  TO  nc(ir) 
ix  =  ix  +  1 

Sigma(ix)  =  SigmaR(ir) 
c(ix)  =  cR(ir) 

Source(ix)  =  SourceR(ir) 
dx(ix)  =  WidthX 

Cel ICtrX ( ix)  =  Xbdy(ir  -  1)  +  (ic  -  .5)  *  WidthX 
NEXT  ic 
NEXT  ir 
ELSE 

SELECT  CASE  DifMeth$ 

CASE  "DD" 

REDIM  Fa(nx,  nk) ,  J(nx,  nk) ,  Jinc(nk) 

CASE  “LD",  “SC”,  “LC",  *  LN " ,  "SA",  “LA",  "EC" 

REDIM  Fa(nx,  nk) ,  Fx(nx,  nk) ,  J(nx,  nk) ,  Jinc(nk) 

END  SELECT 
END  IF 

IF  jinclb  =  0  THEN  'No  incident  current,  so  don't  worry  about  type 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Jinc(k)  =  0 
NEXT  k 
ELSE 

SELECT  CASE  tinclb 

CASE  -1  'Lambertian  incident  current 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Jinc(k)  =  mu(k)  *  4  *  jinclb 
NEXT  k 

CASE  0  'Isotropic  Surface  Source 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Jinc(k)  -  2  *  jinclb 
NEXT  k 

CASE  0  TO  1  'Collimated  Beam 

muO  =  tinclb 
kO  =  nk  +  1 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Jinc ( k)  =0 

IF  muO  <=  mu(k)  THEN  kO  =  k 
NEXT  k 

'kO  is  index  of  smallest  mu  (in  quadrature  set)  <=  muO 
IF  kO  >  nk  THEN 

'muO  >  largest  mu(k) 

Jinc(nk)  =  jinclb  *  (2  /  w(nk) ) 

ELSEIF  muO  =  mu(kO)  THEN 
'muO  =  some  mu(k) 

Jinc(kO)  =  jinclb  *  (2  /  w(k0) ) 
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ELSEIF  kO  =  nk  /  2  +  1  THEN 
'muO  <  smallest  pos .  mu(k) 

Jinc(kO)  =  jinclb  *  (2  /  v;(kO)) 

ELSE 

'mu(kO-l)  <  muO  <  mu(kO)  so  interpolate 
jinc(kO)  =  jinclb  *  (muO  -  mu(kO  -  1) )  /  (mu(kO)  -  mu(kO  -  1)) 
Jinc(kO)  =  Jinc(kO)  *  (2  /  w(kO))  *  (mu(kO)  /  muO) 

J inc ( kO  -  1)  =  jinclb  *  (muO  -  mu(kO))  /  (mu(kO  -  1)  -  mu(kO)) 

J  inc  ( kO  -  1)  =  J  inc  ( kO  -  1)  *  (2  /  w(kO  -  1))  *  (mu(kO  -  1)  /  muO) 
END  IF 
END  SELECT 
END  IF 

IF  jincrb  =  0  THEN  'No  incident  current,  don't  worry  about  type 

FOR  k  =  1  TO  nk  /  2 
Jinc(k)  =  0 
NEXT  k 
ELSE 

SELECT  CASE  tincrb 

CASE  -1  'Lambertian  incident  current 

FOR  k  =  1  TO  nk  /  2 

Jinc(k)  =  mu(k)  *  4  *  jincrb 
NEXT  k 

CASE  0  ' Isotropic  Surface  Source 

FOR  k  =  1  TO  nk  /  2 

Jinc(k)  =  -2  *  jincrb 
NEXT  k 

CASE  0  TO  1  'Collimated  Beam 

muO  =  -tincrb 
kO  =  0 

FOR  k  =  1  TO  nk  /  2 
Jinc(k)  =  0 

IF  muO  >=  mu (k)  THEN  kO  =  k 
NEXT  k 

'kO  is  index  of  least  negative  mu  (in  quadrature  set)  <=  muO 
IF  kO  <  1  THEN 

'muO  <  most  neg .  mu(k) 

Jinc(l)  =  -jincrb  *  (2  /  w(l)) 

ELSEIF  muO  =  mu(kO)  THEN 
'muO  =  some  mu(k) 

Jinc(kO)  =  -jincrb  *  (2  /  w(kO)) 

ELSEIF  kO  =  nk  /  2  THEN 
'muO  >  least  neg.  mu(k) 

Jinc(kO)  =  -jincrb  *  (2  /  w(kO)) 

ELSE 

'mu(kO)  <  muO  <  mu(kO+l)  so  interpolate 
jinc(kO)  =  -jincrb  *  (2  /  w(kO))  *  (mu(kO)  /  muO) 

Jinc(kO)  =  Jinc(kO)  *  (muO  -  mu(kO  +  1))  /  (mu(kO)  -  mu(kO  +  1)) 
Jinc(kO  +  1)  =  -jincrb  *  (2  /  w(kO  +  1))  *  (mu(kO  +  1)  /  muO) 

Jinc(kO  +  l)  =  Jinc(kO  +  1)  *  (muO  -  mu(kO))  /  (mu(kO  +  1)  -  mu(kO)) 
END  IF 
END  SELECT 
END  IF 
RETURN 

SnSol ve : 
iter%  =  0 

GOSUE5  UpdateSourceTermArrays 
DO 

DO 

iter%  =  iter%  +  1 
GOSUB  Enter Left Edge 
GOSUB  WalkRight 
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GOSUB  EnterRightEdge 
GOSUB  WalkLeft 

GOSUB  UpdateF luxe sAndTest Convergence 
GOSUB  UpdateSourceTermArrays 
LOOP  UNTIL  Converged%  OR  (iter%  >=  Itermax%) 

IF  Converged%  THEN 
EXIT  DO 
ELSE 
BEEP 

PRINT  "Not  converged.  MaxChangeObs  =  ";  MaxChangeObs; 

PRINT  "  after";  iter%;  "iterations." 

IF  Choose (" Per form  additional  iterations?  (y/n) :  ")  THEN 
DO 

INPUT  “Additional  number  to  do:  ",  IterAddit ional% 
LOOP  UNTIL  IterAdditional%  >  0 
Itermax%  =  Itermaxl  +  IterAdditional% 

ELSE 

EXIT  DO 
END  IF 
END  IF 
LOOP 
RETURN 

EnterLef tEdge : 

SELECT  CASE  tlb 

CASE  0  'Vacuum 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
J  (  0  ,  k)  =  Jinc (k) 

NEXT  k 

CASE  1  'Symmetry 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 

J(0,  k)  =  Jinc (k)  -  J(0,  nk  -  k  +  1) 

NEXT  k 

CASE  0  TO  1  'Specular  Albedo 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 

J(0,  k)  =  Jinc(k)  -  tlb  *  J(0,  nk  -  k  +  1) 

NEXT  k 

CASE  -1  TO  0  'Grey  Albedo 

Jminus  =  0 
FOR  k  =  1  TO  nk  /  2 

Jminus  =  Jminus  -  J(0,  k)  *  w(k) 

NEXT  k 

Jminus  =  Jminus  *  .5 

Fplus  =  Jminus  *  4  *  ABS(tlb) 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 

J(0,  k)  =  Jinc(k)  +  Fplus  *  mw(k) 

NEXT  k 
END  SELECT 
RETURN 

WalkRight : 

'lprint  "Starting  WalkRight,  iter%  =”;iter% 

SELECT  CASE  DifMeth$ 

CASE  “DD" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

StepDD  J ( ix  -  1,  k) ,  J(ix,  k),  Sa(ix),  dx(ix),  Sigma(ix), 
mu ( k ) ,  Fa ( ix ,  k } 

NEXT  ix 
NEXT  k 
CASE  "LD" 

'  LD  uses  same  Step  algorithm  as  LN,  but  using  a  Pade  expansion 
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'for  exp(-eps)  in  MakeP  converts  implicitly  to  LD 
'Note  this  method  involves  NO  Rotational  Fixup 
FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J(ix  -  1,  k)  /  mu(k) 

StepLN  Sigma(ix),  dx(ix),  mu(k),  Fin,  Sa(ix),  Sx(ix), 
Fa ( ix,  k) ,  Fx ( ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "SC" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J ( ix  -  1,  k)  /  mu(k) 

StepSC  Sigma(ix),  dx(ix),  mu(k),  Fin,  Sa(ix),  Sx(ix), 
Fa ( ix,  k) ,  Fx ( ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "LC" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J ( ix  -  1,  k)  /  mu(k) 

StepLC  Sigma ( ix) ,  dx(ix),  mu(k).  Fin,  Sa(ix),  Sx(ix), 
Fa ( ix,  k) ,  Fx ( ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "LN" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J ( ix  -  1,  k)  /  mu(k) 

StepLN  Sigma(ix),  dx(ix) ,  mu(k) ,  Fin,  Sa(ix),  Sx(ix), 
Fa ( ix,  k) ,  Fx ( ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "SA" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J(ix  -  1,  k)  /  mu(k) 

StepSA  Sigma(ix),  dx(ix),  mu(k).  Fin,  Sa(ix),  Sx(ix), 
Fa  (  ix,  k)  ,  Fx  ( ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "LA" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J(ix  -  1,  k)  /  mu(k) 

StepLA  Sigma(ix),  dx(ix),  mu(k).  Fin,  Sa(ix),  Sx(ix), 
Fa  (  ix,  k)  ,  Fx  (  ix,  k) 

J(ix,  k)  =  Fout  *  mu(k) 

'lprint  j(ix-l,k),  j(ix,k),  Sa(ix) 

'lprint  (sigma(ix)  *  dx(ix)  *  0.5  /  abs(mu(k))),  mu(k) 

Fa ( ix, k) 

NEXT  ix 
NEXT  k 
CASE  "EC" 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
FOR  ix  =  1  TO  nx 

Fin  =  J ( ix  -  1,  k)  /  mu(k) 


Fout , 


Fout , 


Fout , 


Fout , 


Fout , 


Fout , 
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StepEC  Sigma(ix),  dx(ix),  mu(k),  Fin,  Sa(ix),  Sx(ix),  Fout, 

Fa ( ix,  k) ,  Fx ( ix,  k) 

J(ix,  k)  =  Fout  *  rau(k) 

NEXT  ix 
NEXT  k 
END  SELECT 
RETURN 

EnterRightEdge : 

SELECT  CASE  trb 

CASE  0  'Vacuum 

FOR  k  =  1  TO  nk  /  2 
J(nx,  k)  =  Jinc(k) 

NEXT  k 

CASE  1  'Symmetry 

FOR  k  =  1  TO  nk  /  2 

J (nx,  k)  =  Jinc(k)  -  J(nx,  nk  -  k  +  1) 

NEXT  k 

CASE  0  TO  1  'Specular  Albedo 

FOR  k  =  1  TO  nk  /  2 

J (nx,  k)  =  Jinc(k)  -  trb  *  J(nx,  nk  -  k  +  1) 

NEXT  k 

CASE  -1  TO  0  'Grey  Albedo 

Jplus  =  0 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Jplus  =  Jplus  +  J(nx,  k)  *  w(k) 

NEXT  k 

Jplus  =  Jplus  *  .5 

Fminus  =  Jplus  *  4  *  ABS(trb) 

FOR  k  =  1  10  nk  /  2 

J(nx,  k)  =  Jinc(k)  +  Fminus  *  mu(k) 

NEXT  k 
END  SELECT 
RETURN 

WalkLef t : 

'lprint  "Starting  WalkLef t,  iter%  =";iter% 

SELECT  CASE  DifMeth$ 

CASE  "DD" 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 

StepDD  J(ix,  k) ,  J ( ix  -  1,  k) ,  Sa(ix),  dx(ix),  Sigma(ix), 
mu ( k) ,  Fa ( ix,  k ) 

NEXT  ix 
NEXT  k 
CASE  "LD" 

'LD  uses  same  Step  algorithm  as  LN,  but  using  a  Pade  expansion 
'for  exp(-eps)  in  MakeP  converts  implicitly  to  LD 
'Note  this  method  involves  NO  Rotational  Fixup 
FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J ( ix,  k)  /  mu (k) 

StepLN  Sigma(ix),  dx(ix),  -mu(k).  Fin,  Sa(ix),  -Sx(ix),  Fout, 
Fa ( ix,  k) ,  Fx \ ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative  of  mu  and  Sx 
Fx ( ix,  k)  =  -Fx ( ix,  k)  'and  also  get  back  negative  of  Fx 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "SC" 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
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mu  (k) 


-mu(k),  Fin,  Sa(ix),  -Sx(ix),  Fout, 


Fin  =  J(ix,  k)  /  mu(k) 

StepSC  Sigroa(ix),  dx(ix),  -mu(k),  Fin,  Sa(ix),  -Sx(ix),  Fout, 
Fa ( ix,  k) ,  Fx ( ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "LC" 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J(ix,  k)  /  mu(k) 

StepLC  Sigma(ix),  dx(ix),  -mu(k),  Fin,  Sa(ix), 

Fa ( ix,  k)  ,  Fx ( ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 
J(ix  -  1,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  “LN“ 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J(ix,  k)  /  mu(k) 

StepLN  Sigma(ix),  dx(ix),  -mu(k),  Fin,  Sa(ix), 

Fa(ix,  k) ,  Fx(ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  -SA” 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J(ix,  k)  /  mu(k) 

StepSA  Sigma(ix),  dx(ix),  -mu(k),  Fin,  Sa(ix), 

Fa ( ix,  k)  ,  Fx ( ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 


of  mu  and  Sx 


of  mu  and  Sx 
ive  of  Fx 


-mu(k),  Fin,  Sa(ix),  -Sx(ix),  Fout, 


of  mu  and  Sx 
ive  of  Fx 


-mu(k),  Fin,  Sa(ix),  -Sx(ix),  Fout, 


mu  (k) 


Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 

NEXT  ix 
NEXT  k 
CASE  "LA" 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J(ix,  k)  /  mu(k) 

StepLA  Sigma(ix),  dx(ix),  -rou(k).  Fin,  Sa(ix), 
Fa ( ix,  k) ,  Fx( ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 
J ( ix  -  1,  k)  =  Fout  *  mu(k) 

' lprint  j(ix-l,k),  j(ix,k),  Sa(ix) 

'lprint  (sigma(ix)  *  dx(ix)  *  0.5  /  abs(mu(k))) 


of  mu  and  Sx 
ive  of  Fx 


Sa ( ix) 

*  0.5  /  abs (mu ( k) ) ) 


Fa ( ix , k ) 

NEXT  ix 
NEXT  k 
CASE  "EC" 

FOR  k  =  1  TO  nk  /  2 

FOR  ix  =  nx  TO  1  STEP  -1 
Fin  =  J(ix,  k)  /  mu(k) 

StepEC  Sigma(ix),  dx(ix),  -mu(k).  Fin,  Sa(ix), 
Fa ( ix,  k) ,  Fx ( ix,  k) 

'Walk  leftward,  so  reflect  cell,  hence,  use  negative 
Fx(ix,  k)  =  -Fx(ix,  k)  'and  also  get  back  negat 


-Sx ( ix) ,  Fout , 

of  mu  and  Sx 
ive  of  Fx 


,  mu ( k ) , 


-Sx ( ix) ,  Fout , 

of  mu  and  Sx 
ive  of  Fx 


Fout  *  mu(k) 


J ( ix  -  1 ,  k)  = 
NEXT  ix 
NEXT  k 
END  SELECT 
RETURN 


UpdateFluxesAndTestConvergence : 
MaxChangeObs  =  0 
SELECT  CASE  DifMeth$ 

CASE  "DD" 

Sum  =  0 

FOR  k  =  1  TO  nk 

Sum  =  Sum  +  w(k)  *  J(0,  k)  /  mu(k) 
NEXT  k 

BdyFlux  =  Sum  /  2 
FOR  ix  =  1  TO  nx 
Sum  =  0 

FOR  k  =  1  TO  nk 

Sum  =  Sum  +  w{k)  *  Fa(ix,  k) 
NEXT  k 


FluxNewA  =  Sum  /  2 

IF  (ABS (FluxNewA)  <=  Change)  AND  ( ABS ( FluxA ( ix) )  <=  Change) 
ChangeObs  =  MAX (ABS (FluxNewA) ,  ABS (FluxA( ix) ) ) 

ELSE 


ChangeObs  =  ABS (FluxNewA  -  FluxA (ix))  *  2  / 
FluxA ( ix) ) 

END  IF 

MaxChangeObs  =  MAX (MaxChangeObs ,  ChangeObs) 
FluxA (ix)  =  FluxNewA 
Sum  =  0 

FOR  k  =  1  TO  nk 

Sum  =  Sum  +  w(k)  *  J(ix,  k)  /  mu(k) 

NEXT  k 

FluxX(ix)  =  (Sum  /  2)  -  BdyFlux 
BdyFlux  =  Sum  /  2 
NEXT  ix 
CASE  "LN" 

FOR  ix  =  1  TO  nx 
SumA  =  0 
SumX  =  0 


FOR  k  =  1  TO  nk 

SumA  =  SumA  +  w(k)  *  Fa(ix,  k) 
SumX  =  SumX  +  w(k)  *  Fx(ix,  k) 
NEXT  k 


(FluxNewA  + 


FluxNewA  =  SumA  /  2 
FluxNewX  =  SumX  /  2 

IF  (ABS (FluxNewA)  <=  Change)  AND  ( ABS ( FluxA ( ix) )  <=  Change) 
ChangeObs  =  MAX (ABS ( FluxNewA)  ,  ABS ( FluxA ( ix) ) ) 

ELSE 

ChangeObs  =  ABS(FluxNewA  -  FluxA(ix))  *  2  /  (FluxNewA  + 
FluxA ( ix) ) 

END  IF 

MaxChangeObs  =  MAX (MaxChangeObs ,  ChangeObs) 

FluxA (ix)  =  FluxNewA 
FluxX(ix)  =  FluxNewX 

IF  fixup%  AND  ABS ( FluxX ( ix) )  >  FluxA(ix)  THEN 
FluxX(ix)  =  SGN (FluxX ( ix) )  *  FluxA(ix) 


END  IF 
NEXT  ix 

CASE  " LD“ ,  "SC”,  "LC",  "SA",  "LA",  "EC" 

FOR  ix  =  1  TO  nx 


SumA  =  0 


THEN 


THEN 
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SumX  =  0 

FOR  k  =  1  TO  nk 

SumA  =  SumA  +  w(k)  *  Fa(ix,  k) 

SumX  =  SumX  +  w(k)  *  Fx(ix,  k) 

NEXT  k 

FluxNewA  =  SumA  /  2 
FluxNewX  =  SumX  /  2 

IF  (ABS ( FluxNewA)  <=  Change)  AND  (ABS ( FluxA ( ix) )  <=  Change)  THEN 
ChangeObs  =  MAX (ABS ( FluxNewA) ,  ABS (FluxA ( ix) ) ) 

ELSE 

ChangeObs  =  ABS(FluxNewA  -  FluxA(ix))  *  2  /  (FluxNewA  + 

FluxA ( ix) ) 

END  IF 

MaxChangeObs  =  MAX (MaxChangeObs ,  ChangeObs) 

FluxA (ix)  =  FluxNewA 
FluxX(ix)  =  FluxNewX 
NEXT  ix 
END  SELECT 
LOCATE  15,  1 

PRINT  “  Working  *;  DifMeth$;  TAB(15); 
timiters  =  (TIMER  -  StartClock)  /  iter%  /  60 
PRINT  USING  ####.####" ;  timiters; 

PRINT  "  min/ iter  "; 

IF  ( iter%  =  1)  AND  (iterl%  >  1)  THEN 
LOCATE  16,  15 

timiters  =  timiters  *  iterl% 

PRINT  USING  "Est  ####.####";  timiters; 

PRINT  “  min  to  go  “ 

END  IF 
LOCATE  18,  1 

PRINT  "  After  Iteration";  iterl%;  "  MaxChangeObs  =  ";  MaxChangeObs 1 ; 

M 

PRINT  “  “;  iter%;  *  =  ";  MaxChangeObs ;  " 

iterl%  =  iter% 

MaxChangeObs 1  =  MaxChangeObs 
PRINT 

Converged%  =  (MaxChangeObs  <  Change) 

RETURN 


UpdateSourceTermArrays : 

SELECT  CASE  DifMeth$ 

CASE  "DD" 

FOR  ix  =  1  TO  nx 

Sa(ix)  =  c(ix)  *  Sigma (ix)  *  FluxA(ix)  +  Source (ix) 


NEXT  ix 
CASE  “LC" 

FOR  ix  =  1  TO  nx 

Sa(ix)  =  c(ix)  *  Sigma(ix)  *  FluxA(ix)  +  Source(ix) 
Sx(ix)  =  c(ix)  *  Sigma(ix)  *  FluxX(ix) 

'Implement  Negative  Source  Fixup  as  required  for  LC 
IF  ABS ( Sx ( ix ) )  >  Sa(ix)  THEN 
Sx(ix)  =  SGN(Sx(ix))  *  Sa(ix) 


END  IF 
NEXT  ix 

CASE  "LD",  “LN",  "SC",  "SA",  "LA",  "EC" 

FOR  ix  =  1  TO  nx 

Sa(ix)  =  c(ix)  *  Sigma(ix)  *  FluxA(ix)  +  Source(ix) 
Sx(ix)  =  c(ix)  *  Sigma(ix)  *  FluxX(ix) 

NEXT  ix 
END  SELECT 
RETURN 
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CalculateResults : 

REDIM  JPlusBdy (nr ) ,  JMinusBdy (nr ) ,  JNetBdy(nr),  FluxBdy(nr) 

REDIM  FluxAveR (nr ) ,  FluxXR(nr) 
ix  =  0 
ir  =  0 

GOSUB  CalculateRegionBoundaryResults 
'left  bdy  of  first  region 
FOR  ir  =  1  TO  nr 

GOSUB  CalculateRegionAverageScalarFlux 
'ir'th  region 

GOSUB  CalculateRegionBoundaryResults 
'right  bdy  of  ir'th  region 
NEXT  ir 
RETURN 

CalculateRegionBoundaryResults : 

Sum  =  0 

FOR  k  =  nk  TO  nk  /  2  +  1  STEP  -1 
Sum  =  Sum  +w(k)  *J(ix,  k) 

NEXT  k 

JPlusBdy (ir)  =  Sum  /  2 
Sum  =  0 

FOR  k  =  1  TO  nk  /  2 

Sum  =  Sum  +w(k)  *  J(ix,  k) 

NEXT  k 

JMinusBdy ( ir)  =  -Sum  /  2 

JNetBdy(ir)  =  JPlusBdy(ir)  -  JMinusBdy ( ir ) 

Sum  =  0 

FOR  k  =  1  TO  nk 

Sum  =  Sum  +  w(k)  *  J(ix,  k)  /  mu(k) 

NEXT  k 

FluxBdy(ir)  =  Sum  /  2 
RETURN 

CalculateRegionAverageScalarFlux: 

SumA  =  0 
SumX  =  0 
Suml  =  0 

FOR  ic  =  1  TO  nc(ir) 
ix  =  ix  +  1 

SumA  =  SumA  +  FluxA(ix) 

SumX  =  SumX  +  FluxX(ix) 

Suml  =  Suml  +  (2  *  ic  -  1)  *  FluxA(ix) 

NEXT  ic 

FluxAveR (ir)  =  SumA  /  nc(ir) 

FluxXR(ir)  =  -3  *  SumA  /  nc(ir)  +  (SumX  +  3  *  Suml)  /  nc(ir)  *  2 
RETURN 

VerifyRegionsByBalanceEquation : 

MaxRegrerr  =  0 ! 

RegBalFlag$  =  “ " 

Regflag%  =  False 
FOR  ir  =  1  TO  nr 

SumBL  =  ( JNetBdy ( ir )  -  JNetBdy(ir  -  1))  /  (Xbdy(ir)  -  Xbdy(ir  - 
+  SigmaR(ir)  *  FluxAveR (ir) 

SumBR  =  (cR(ir)  *  SigmaR(ir)  *  FluxAveR(ir)  +  SourceR(ir)) 

IF  (SumBL  <>  0)  OR  (SumBR  <>  0)  THEN 

Regrerr  =  ABS (SumBL  -  SumBR)  *  2  /  (ABS (SumBL)  +  ABS (SumBR) ) 
ELSE 

Regrerr  =  0 
END  IF 

IF  Regrerr  <  2  THEN  MaxRegrerr  =  MAX (MaxRegrerr ,  Regrerr) 


IF  MaxRegrerr  >  Change  THEN 
Regflag%  =  True 

RegBalFlag$  =  LTRIM$ (RegBalFlag$  +  STR$(ir)) 

END  IF 
NEXT  ir 
RETURN 

LPr  intRegionSummary : 

SELECT  CASE  DevEcho$ 

CASE  " LPT1 “ 

OPEN  " LPT1 : “  FOR  OUTPUT  AS  #4 
CASE  "SCREEN" 

OPEN  " SCRN : "  FOR  OUTPUT  AS  #4 
CASE  “OUT" 

OPEN  Outf ile$  FOR  APPEND  AS  #4 
END  SELECT 
SELECT  CASE  DifMeth$ 

CASE  "DD“ 

IF  fixup%  THEN 

PRINT  #4,  "DDF  Negative  Flux  Fixup  Enabled" 

ELSE 

PRINT  #4,  “DD  No  Flux  Fixup" 

END  IF 
CASE  “LC" 

PRINT  #4,  "LC  Source  Rotation  (Sx<=Sa)  Fixup  ALWAYS  Enabled" 

CASE  "LN" 

IF  f ixup%  THEN 

PRINT  #4,  “LNF  Scalar  Flux  Rotation  (PhiX<PhiA)  Fixup  Enabled" 
ELSE 

PRINT  #4,  "LN  No  Flux  Fixup" 

END  IF 
CASE  "EC" 

IF  SwitchEC%  THEN 

PRINT  #4,  "EC  to  LC  Enabled  (beta  <="; 

ELSE 

PRINT  #4,  "EC  Normal  (beta  >  "; 

END  IF 

PRINT  #4,  USING  “  ##  .  «  ;  SwECtoLCSet Pt ; 

PRINT  #4,  * ) " 

CASE  ELSE 

PRINT  #4,  DifMeth$ 

END  SELECT 

PRINT  *4,  "Execution  Time  :  "; 

PRINT  #4,  USING  "###.####";  ExecMin; 

PRINT  #4,  "  min  " 

IF  Converged%  THEN  PRINT  #4,  “Converged  "; 

PRINT  #4,  "After";  iter%;  "iterations,  MaxChangeObs  =";  MaxChangeObs 
PRINT  #4, 

PRINT  #4,  "  J  plus",  "  J  minus",  "  J  net",  "  Boundary  Flux" 

PRINT  #4,  "Region  #",  "  flux  ave",  "  flux  x-moment" 

PRINT  #4, 
ir  =  0 

PRINT  #4,  “ " , 

PRINT  #4,  USING  "  +#  .  #####*~/VA/'  *;  JPlusBdy  (  ir )  ;  JMinusBdy ( ir ) ; 

PRINT  #4,  USING  *  +# .  ######----  JNetBdy(ir);  FluxBdy(ir) 

FOR  ir  =  1  TO  nr 
PRINT  #4,  ir, 

PRINT  #4,  USING  "  +#  .  ■;  FluxAveR ( ir ) ;  FluxXR(ir) 

PRINT  #4,  "", 

PRINT  #4,  USING  "  +  #.######  AA/VA  *;  JPlusBdy  (  i  r  )  ;  JMinusBdy ( ir ) ; 
PRINT  #4 ,  USING  ’  +  #.######A~^  ";  JNetBdy(ir);  FluxBdy(ir) 

NEXT  ir 
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IF  (MaxMomOrerr  >  MomTol)  THEN 
MomOIntF lag$  =  “  Oth:  VIOLATION" 

END  IF 

IF  (MaxMomlrerr  >  MomTol)  THEN 
MomlIntFlag$  =  “  1st:  VIOLATION" 

END  IF 
PRINT  #4, 

PRINT  #4,  "Moment  Balances:  " ;  MomOIntFlag$ ; 

PRINT  #4,  USING  "  MaxRelErr:  ##.####----■;  MaxMomOrerr; 

PRINT  #4,  “  <  *; 

PRINT  #4,  USING  *  ##  .  ####/'A/'/'  ■  ;  MomTol 
SELECT  CASE  DifMeth$ 

CASE  “DD",  "SC" 

CASE  ELSE 

PRINT  #4,  TAB (18) ;  Moml IntFlag$ ; 

PRINT  #4,  USING  "  ##  .  ####A/V/V"  “  ;  MaxMomlrerr; 

PRINT  #4,  "  <  ■; 

PRINT  #4,  USING  "  ##  .  ####AAAA  •  ;  MomTol 
END  SELECT 

PRINT  #4,  "  Region  Balances:  "; 

IF  Regf lag%  =  True  THEN 

PRINT  #4,  TAB (24) ;  LTRIM$ ( "VIOLATED  in  “  +  RegBalFlag$  +  “:  "); 

ELSE 

PRINT  #4,  "CONSERVED";  TAB(48); 

END  IF 

PRINT  #4,  USING  " ##  .  ####AA AA “ ;  MaxRegrerr; 

PRINT  #4,  "  <  "; 

PRINT  #4,  USING  *  ##  .  ####AAAA "  ;  Change 
PRINT  #4,  :  PRINT  #4, 

CLOSE  #4 
RETURN 

RunSummary : 

SELECT  CASE  DevEcho$ 

CASE  “ LPT1 " 

OPEN  " LPT1 : "  FOR  OUTPUT  AS  #4 
CASE  "SCREEN" 

OPEN  " SCRN : "  FOR  OUTPUT  AS  #4 
CASE  "OUT" 

OPEN  Out  f i le$  FOR  APPEND  AS  #4 
END  SELECT 

MinCh’ngeReg  -  1E+15 
ilotogl  =  -2:  ihitogl  =  0 
FOR  it  =  1  TO  itim 

GOSUB  ToggleDi fMeth 

IF  DifMeth$  =  "EC"  THEN  iec  =  it 

MmChangeReg  =  MIN  (MinChangeReg ,  RErranal  (  it )  ) 

RelRt im ( it )  =  (1000!  *  Rtim(it)) 

RelRErranal ( it )  =  (1000!  *  RErranal ( it ) ) 

NEXT  it 

IF  iec  <  1  OR  iec  >  itim  THEN  iec  =  1 

PRINT  #4,  -************************************************************ 

*************** 

PRINT  #4, 

PRINT  #4,  "  RunSummary:  Relative  Processing  times  and  Region 

MaxChanges " 

PRINT  #4,  “  Method  Exec  Min  RelTime  RegMaxErr 

RelRegErr " 

PRINT  #4, 

ilotogl  =  -2:  ihitogl  =  0 
FOR  it  =  1  TO  itim 

GOSUB  ToggleDi fMeth 
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IF  (DifMeth$  =  *DD"  OR  DifMeth$  =  *LN“)  AND  (fixup%)  THEN 
f ixa$  =  “F* 

ELSE  f ixa$  =  *' 

END  IF 

PRINT  #4,  TAB ( 8 )  ;  LTRIM$ (DifMeth$  +  fixa$);  TAB(14); 

PRINT  #4,  USING  “####.####“;  Rtim(it); 

PRINT  #4,  TAB (26)  ; 

PRINT  #4,  USING  *###.####“;  CSNG (RelRt im ( it )  /  RelRt im ( iec ))  ; 
PRINT  #4,  TAB (37) ; 

PRINT  #4,  USING  "##.### * ;  RErranal ( it ) ; 

PRINT  #4,  TAB (51) ; 

IF  MinChangeReg  <>  0  THEN 

PRINT  #4,  USING  ■###.####«;  CSNG ( Re IRErranal ( i t )  /  (1000  * 
MinChangeReg) ) 

ELSE  PRINT  #4,  * - * 

END  IF 
NEXT  it 

PRINT  #4,  :  PRINT  #4, 

CLOSE  #4 
RETURN 


CreateTKDataFi le : 

SELECT  CASE  CmdProf$ 

CASE  "S" 

DO 

INPUT  "Output  Data  to  a  TK  File?  (  Y/N  ):  ”,  Qtkf$ 

Qtkf$  =  UCASE$ (LEFTS (LTRIM$ (Qtkf $ ) ,  1)) 

LOOP  UNTIL  INSTR ( “ YN" ,  Qtkf$)  AND  LEN (Qtkf $ )  =  1 
CASE  "A" 

END  SELECT 

IF  UCASE$ ( Qtkf  $ )  =  “Y"  THEN 
SELECT  CASE  CmdProf$ 

CASE  "S" 

DO 

INPUT  "TK  File  name?  :  *,  TKfile$ 

LOOP  WHILE  TKfileS  <> 

OPEN  TKf i le$  FOR  APPEND  AS  #3 
CASE  “ A“ 

OPEN  TKf i le$  FOR  APPEND  AS  #3 
END  SELECT 

IF  (DifMethS  =  -DD"  OR  DifMeth$  =  " LN" )  AND  (fixup%)  THEN 
f ixa$  =  "F" 

ELSE  f ixa$  =  "" 

END  IF 

CALL  PrtTK (LTRIM$ ( “cR_"  +  DifMethS  +  fixa$),  cR(),  1,  nr) 

CALL  PrtTK (LTRIMS ( "SigR_"  *  DifMeth$  +  fixa$),  SigmaRO,  1,  nr) 

CALL  PrtTK  (LTRIM$  ( “Src_"  +  DifMethS  +  fixa$),  SourceRO,  1,  nr) 
REDIM  cellsnc(l  TO  nr) 

FOR  ir  =  1  TO  nr 

cellsnc(ir)  =  nc(ir) 

NEXT  ir 

CALL  PrtTK  ( LTRIMS  (  ”nc_*  +  DifMethS  +  fixa$),  cellsncO,  1,  nr) 

CALL  PrtTK ( LTRIMS ( "Xb_“  +  DifMethS  +  fixa$),  Xbdy ( ) ,  0,  nr) 

CALL  PrtTK (LTRIMS ( "Jposb_"  +  DifMethS  +  fixa$),  JPlusBdy ( ) ,  0,  nr) 

CALL  PrtTK (LTRIMS ( "Jnegb_“  +  DifMethS  +  fixaS),  JMinusBdy ( ) ,  0,  nr) 

CALL  PrtTK  (LTRIMS  (  *Jnetb_"  +  DifMethS  +  fixaS),  JNetBdyO,  0,  nr) 

CALL  PrtTK  ( LTRIMS  ("Phib_“  +  DifMethS  4  fixa$),  FluxBdyO,  0,  nr) 
REDIM  Xctr(nr) 

FOR  ir  =  1  TO  nr 

Xctr(ir)  =  Xbdvfir  -  1)  4  (Xbdyfir)  -  Xbdy(ir  -  1))  /  2! 

NEXT  ir 

CALL  PrtTK (LTRIMS ( “Xctr_“  4  DifMethS  4  fixa$),  Xctr(),  1,  nr) 
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CALL  PrtTK  (LTRIM$  ( *  PhiA_"  +  DifMeth?  +  fixa$),  FluxAveRO,  1,  nr) 
CALL  PrtTK (LTRIM$ ( "PhiX_"  +  DifMeth$  +  fixa$),  FluxXRO,  1,  nr) 
'establish  a  global  flux  variable 
ng  =  2  *  nr 

REDIM  XGlbl ( 0  TO  ng),  FlxGlbl(0  TO  ng) 
ig  =  0 :  ix  =  0 :  iy  =  1 
DO 

XGlbl (ig)  =  Xbdy(ix) 

FlxGlbl ( ig)  =  FluxBdy(ix) 
ix  =  ix  +  1 
ig  =  ig  +  1 

IF  ix  <>  (nr  +  1)  THEN 
XGlbl (ig)  =  Xctr(iy) 

FlxGlbl(ig)  =  FluxAveR(iy) 
iy  =  iy  +  1 
ig  =  ig  +  1 
END  IF 

LOOP  UNTIL  ix  =  (nr  +  1) 

CALL  PrtTK (LTRIM$ ( "XG_"  +  DifMeth$  +  fixa$),  XGlbl () ,  0,  ng) 

CALL  PrtTK (LTRIM$ ( "PhiG_"  +  DifMeth?  +  fixa$),  FlxGlbl () ,  0,  ng) 

CALL  PrtTK(LTRIM$ ( "Time_*  +  DifMeth$  +  fixa$),  Rtim(),  itim,  itim) 
PRINT  *3, 

CLOSE  #3 
END  IF 
RETURN 

SUB  AltMenu  STATIC 

SHARED  Problemf ile$,  Outfilo$,  Qoutf$,  TKfile$,  Qtkf$,  DifMeth$, 

TotDi f Meth$ 

SHARED  ilotogl,  ihitogl,  QuadMeth$,  fixup%,  Flfix$,  Itermaxl,  Change 
SHARED  MomTol,  nk,  Ident$,  again%,  NewProblem%,  ResetAl lFluxes$ 

SHARED  SwECtoLCSetPt 
Variables : 

CLS 

PRINT  ****************************************************************** 
***★*★**★•• 

PRINT  "Problem  File  :  * ;  UCASE$ (Problemf ile$) ;  TAB(40);  "nk  Angular 

Ordinates  : “ ;  nk 

PRINT  "Output  file  :  ";  UCASE$ (Outf ile$) ;  TAB(40);  "Quadrature 

Method  :  UCASE$ (QuadMeth$ )  +  LTRIM$ ( STR$ ( nk) ) 

PRINT  "Tk  file  :  ";  UCASE$ (TKf ile$ ) ;  TAB(40);  "Negative  Flux 

Fixups  :  "?  UCASE$ (Flf ix$) 

PRINT  "Max  Iterations  Itermax%;  TAB(40);  "Reset  All  Old  Fluxes  : 

ResetAl lFluxes$ 


PRINT  "Moment  Tolerance  : ” ; 

PRINT  USING  MomTol; 

PRINT  TAB (40) ;  "Solution  Tolerance 
PRINT  USING  "##  .  ;  Change 

PRINT 

PRINT  "  Spatial  Quadrature  Methods  TotDifMeth$ 

PRINT  "  File  Identifier  :  " ;  Ident$ 

PRINT  -***************************************************************** 
PRINT 

PRINT  "Type  to  change  •;  TAB(40);  "Type  to  change  " 

PRINT  " -  - *;  TAB  (40)  ;  " - 


IF  NOT  again%  THEN 

PRINT  "  PF  Problem  File  TAB(40);  “  NK  #  Angular  Ordinates  " 

PRINT  "  O  Data  Output  File";  TAB(40);  "  QM  Quadrature  Method" 

PRINT  "  TK  TK  Solver  Output  File";  TAB(40);  "  FF  Neg  Flux  Fixups 
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ELSE 

PRINT 

PRINT 

PRINT 


Problem  File  TAB (40);  "  NK  #  Angular  Ordinates  “ 

Data  Output  File”;  TAB(40);  “  QM  Quadrature  Method” 

TK  Solver  Output  File";  TAB(40);  "  FF  Neg  Flux  Fixups 


END  IF 

PRINT  “  MI  Max  Iterations  Allowed  TAB(40);  “  TL  Change 
Tolerances " 

PRINT  “  AS  Add  Spatial  Quad  Method  ";  TAB(40);  "  NF  Reset  Flux 

Usage  Toggle" 

PRINT  "  CS  Clear  Spatial  Quad  Methods  ";  TAB(40);  "  FL  DOS  File 

Directory" 

PRINT  "  N  None  of  the  above" 

PRINT 

INPUT  Response$ 

LOCATE  12,  1 
FOR  i  =  1  TO  12 

PRINT  TAB (79 ) ;  '  " 

NEXT  i 

LOCATE  12,  1 

IF  UCASES  (Responses)  =  "FL"  AND  NOT  agam%  THEN 
INPUT  “Path  or  file:  ";  fileinqS 
FILES  (fileinqS) 

PRINT  :  PRINT  "Holding  5  Seconds...  PRINT 
SLEEP  (5) 

ELSEIF  UCASES (Responses)  =  "PF"  AND  NOT  again%  THEN 
INPUT  "Problem  File  :  ";  Problemfile$ 

ELSEIF  UCASES (Responses)  =  *0“  AND  NOT  again%  THEN 
INPUT  “Output  file  (OFF  for  none)  :  ";  OutfileS 
IF  UCASES (OutfileS)  =  "OFF"  THEN 
Qoutf $  =  “N" 

ELSE 

QoutfS  =  “ Y “ 

QtkfS  =  "Y" 
ippos  =  0 

ippos  =  INSTR (Out  file$,  ".")  -  1 

TKfileS  =  LTRIMS (MIDS (OutfileS,  1,  ippos)  +  " . TKD" ) 

END  IF 

ELSEIF  UCASES (Responses)  =  "TK"  AND  NOT  again%  THEN 
INPUT  "Tk  file  (OFF  for  none)  :  ";  TKfileS 
QtkfS  =  “Y" 

IF  UCASES (TKfileS)  =  "OFF"  THEN 
QtkfS  =  "N" 

END  IF 

ELSEIF  UCASES (Responses)  =  "MI"  THEN 
DO 

INPUT  "Maximum  Allowable  Iterations  ";  Itermax% 

LOOP  UNTIL  Itermax%  >  1 
ELSEIF  UCASES (Responses)  =  "NK"  THEN 
SELECT  CASE  QuadMethS 
CASE  "S" 

DO 

PRINT  "For  Single-Range  Gauss  Sn,  n  is  total  #  of  mu's." 

PRINT  "Supported  Orders  are  n  =  2,  4,  6,  8,  10  and  12." 

PRINT 

INPUT  "Enter:  n  =  ",  nk 

LOOP  UNTIL  (nk  >  =  2)  AND  (nk  <=  12)  AND  (nk  MOD  2=0) 

CASE  "D" 

DO 

PRINT  "For  Double-Range  Gauss  Sn ,  n  is  #  of  mu's  in  each  range." 

PRINT  "Supported  Orders  are  n  =  1,  2,  3,  4,  6,  8,  10  and  Xl." 

PRINT 
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INPUT  “Enter:  N  =  " ,  nk 

LOOP  UNTIL  nk  >=  1  AND  nk  <=  12  AND  ( (nk  MOD  2=0)  OR  nk  <  4) 

CASE  “M" 

DO 

PRINT  “For  Composite  Midpoint  Rule  Sn,  n  is  total  #  of  mu's. “ 
PRINT  "The  mu's  are  at  the  center  of  equally  sized  intervals," 
PRINT  “  and  are  given  equal  weights." 

PRINT  "Supported  Orders  are  n  =  2,  4,  6,  ..." 

PRINT 

INPUT  "Enter:  N  =  ",  nk 
LOOP  UNTIL  nk  >  0  AND  ( (nk  MOD  2)  =  0) 

END  SELECT 

ELSEIF  UCASE$ (Response$ )  =  “QM*  THEN 
DO 

PRINT  "Select  Angular  Quadrature  Type:":  PRINT 
PRINT  “  S  -  Single  Range  Gauss-Legendre" 

PRINT  "  D  -  Double  Range  Gauss-Legendre" 

PRINT  "  M  -  Composite  Midpoint  Rule":  PRINT 
INPUT  “Choice?  (  S,D,M  ):  ",  QuadMeth$ 

QuadMethS  =  UCASE$ (LEFTS (LTRIM$ (QuadMeth$)  ,  1)) 

LOOP  UNTIL  INSTR ( " SDM" ,  QuadMethS)  AND  LEN (QuadMethS )  =  1 
ELSEIF  UCASES (Responses)  =  “FF"  THEN 
DO 

INPUT  "Provide  Negative  Flux  Fixups  (Enter  Y  or  N  only)  FlfixS 
LOOP  UNTIL  (UCASES (FlfixS)  =  "Y")  OR  (UCASES ( Fl f ix$ )  =  "N" ) 

IF  UCASES (FlfixS)  =  "Y"  THEN 

FlfixS  =  "YES":  f ixup%  =  True 
ELSE 

FlfixS  =  "NO":  fixup%  =  False 
END  IF 

ELSEIF  UCASES (Responses)  =  “TL"  THEN 
DO 

INPUT  "Tolerance  for  Solution  Convergence  (0  to  1)  :  *;  Change 

LOOP  UNTIL  Change  >0# 

PRINT 

DO 

INPUT  "Tolerance  for  Moment  Balance  Comparisons  (0  to  1)  :  MomTol 

LOOP  UNTIL  Change  >0# 

ELSEIF  UCASES (Responses)  =  "NF"  THEN 
DO 

INPUT  "Reset  all  Fluxes  ?  (  Y/N  ):  ",  a$ 
aS  =  UCASES (LEFTS (LTRIMS (a$) ,  1)) 

LOOP  UNTIL  INSTR ( " YN"  ,  a$)  AND  LEN (a$ )  =  1 
IF  a$  =  "Y"  THEN 

Res e tAl IF luxes $  =  "YES" 

ELSE 

Reset A1 lFluxesS  =  "NO” 

END  IF 

ELSEIF  UCASES (Responses)  =  "AS"  THEN 
IF  TotDi fMeth$  =  "  NONE"  THEN  TotDifMethS  =  ■  " 

DO 

PRINT  "  Add  Spatial  Quadrature  Method(s)  :" 

PRINT 

PRINT  "  DD  -  Diamond  Difference  MB  -  DD,  LD, LC , SA, LA, and 

EC" 

PRINT  *  LD  -  Linear  Discontinuous  EL  -  Entire  List  " 

PRINT  “  SC  -  Step  Characteristic" 

PRINT  "  LC  -  Linear  Characteristic" 

PRINT  "  LN  -  Linear  Nodal" 

PRINT  "  SA  -  Step  Adaptive" 

PRINT  "  LA  -  Linear  Adaptive" 

PRINT  "  EC  -  Exponential  Characteristic" 
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PRINT 

INPUT  “Choice?  (  DD, LD, SC, LC, LN, SA, LA, EC, MB,  EL  ):  “,  DifMeth$ 
DifMethS  =  UCASE$ (LEFT$ (LTRIM$ (DifMeth$ )  ,  2)) 

LOOP  UNTIL  INSTR ( " ELMBDDLDSCLCLNSALAEC “ ,  DifMethS)  AND  LEN(DifMeth$) 

2 

SELECT  CASE  DifMeth$ 

CASE  “MB" 

TotDi fMeth$  =  "  DD  LD  LC  SA  LA  EC  “ 

CASE  "EL" 

TotDi fMeth$  =  “  DD  LD  SC  LC  LN  SA  LA  EC  “ 

CASE  ELSE 

IF  INSTR (TotDi fMeth$ ,  DifMeth$)  THEN 
ELSE 

TotDifMeth$  =  TotDifMeth$  +  LTRIM$ (DifMeth$  +  "  ") 

END  IF 
END  SELECT 

ELSEIF  UCASE$ (Responses)  =  “CS“  THEN 
TotDi fMeth$  =  *  NONE" 

ELSETF  UCASE$ (Responses)  =  “N“  AND  TotDifMethS  <>  '  NONE"  THEN 
IF  INSTR (TotDifMethS,  “EC")  THEN 
CALL  SetSwitchECtoLC 
LOCATE  12,  1 
FOR  i  =  1  TO  12 

PRINT  TAB (79) ;  “  " 

NEXT  i 

LOCATE  12,  1 

END  IF 

IF  Problemf ile$  =  ""  THEN 

LOCATE  12,  1:  PRINT  “WARNING:  No  Problem  File  is  Designated" 

SLEEP  (3) 

GOTO  Variables 
END  IF 

IF  UCASES (Outf ile$)  =  “OFF"  THEN 

LOCATE  12,  1:  PRINT  “WARNING:  No  Output  File  is  Designated” 

PRINT 

DO 

INPUT  "Proceed  ?  (  Y/N  ):  ",  a$ 
a$  =  UCASES  (LEFTS  (LTRIMS  (a$)  ,  D) 

LOOP  UNTIL  INSTR (“YN",  a$)  AND  LEN (a$ )  =  1 
LOCATE  12,  1 
FOR  i  =  1  TO  10 

PRINT  TAB (79) ;  “  " 

NEXT  i 

LOCATE  12,  1 

IF  a$  =  "N"  THEN  GOTO  Variables 
END  IF 
GOTO  Bottom 
END  IF 

GOTO  Variables 
Bottom : 

'  Open  files  as  Necessary 
IF  NOT  again%  THEN 
OPEN  Problemf ile$  FOR  INPUT  AS  #1 
'check  for  correct  file: 

INPUT  #1,  IdentS 

LOCATE  9,  38:  PRINT  IdentS:  LOCATE  12,  1 
DO 

INPUT  "Correct  File?  (  Y/N/Quit  ):  “,  a$ 
a$  =  UCASES  (LEFTS  (LTRIMS  (a$)  ,  D) 

LOOP  UNTIL  INSTR ( ’ YNQ" ,  a$)  AND  LEN (a$ )  =  1 
SELECT  CASE  a$ 

CASE  "N" 
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CLOSE  #1 
Problemfile?  = 

Ident  $  =  "" 

GOTO  Variables 
CASE  "Q" 

END 

CASE  “Y" 

CLS 

CASE  ELSE 
BEEP 

PRINT  “ERROR:  Unsupported  choice  in  OpenProblemFi le . “ 

STOP 
END  SELECT 
END  IF 

'  Assign  First  String  toggles  and  Corresponding  Diff  Method 
ilotogl  =  1:  ihitogl  =  3 

DifMeth$  =  LTRIM$ (MID$ (TotDifMeth$ ,  ilocogl,  ihitogl)) 

DifMeth$  =  UCASE$ (LEFTS (LTRIM$ (DifMeth$)  ,  2)) 

END  SUB 

SUB  BeginSets 

SHARED  Problemf i le$ ,  DifMeth$,  QuadMeth$,  fixup%,  Flfix$,  Itermax% 
SHARED  Change,  MomTol,  nk,  TotDifMeth$,  Outfile$,  Qoutf?,  TKfile$ 
SHARED  Qtkf$,  Ident$,  again%,  NewProblem%,  ResetA1 lFluxes? 

SHARED  MomOIntFlag?,  MaxMomOrerr,  MomlIntFlag$,  MaxMomlrerr 
SHARED  SwitchEC% ,  SwECtcLCSet Pt 

’  This  provides  the  user  with  input  values 
'  for  an  initial  run. 

t 

Problemfile$  =  “ a : test  1 . IPT“ 

DifMeth$  =  "DD" 

TotDifMeth$  =  “  NONE" 

QuadMeth$  =  "S" 
fixup%  =  False 
Flf ix$  =  "NO" 

Itermax%  =  150 
Change  =  .00001# 

MomTol  =  IE-10 
nk  =  8 
Ident$  =  -  " 
again%  =  False 
NewProblem%  =  True 
SwitchEC%  =  False 
SwECtoLCSet  Pt  =  .0005 
ResetAl lFluxes$  =  "YES" 

Outfile$  =  "off" 

Qoutf?  =  "n" 

TKf i le$  =  "off" 

Qtkf $  =  "n" 

Mom0IntFlag$  =  "  0th:  CONSERVED":  MaxMomOrerr  =  0 
Moml IntFlaa$  =  "  1st:  CONSERVED":  MaxMomlrerr  =  0 
END  SUB 

DEFINT  A-H,  J,  L-M,  O-Z 
FUNCTION  Choose  (prompt?)  STATIC 
CONST  False  =  0 
CONST  True  =  NOT  False 
DO 

PRINT  prompt?; 

INPUT  * " ,  ch$ 

ch$  =  UCASE? (LEFT? (LTRIM? (ch$)  ,  1)) 
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LOOP  UNTIL  INSTR ( “YN10 " ,  ch$)  AND  LEN ( ch$ )  = 

IF  ch$  =  “ Y “  OR  ch$  =  "1“  THEN 
Choose  =  True 
ELSE 

Choose  =  False 
END  IF 

END  FUNCTION 
DEFSNG  H 

DEFDBL  A-G,  J,  L-M,  O-Z 
FUNCTION  DdxG  (x)  STATIC 
'  d/dx  g(x)  for  Newton's  Method  Root  Solving  beta  in  EC 
DdxG  =  1  /  x  2  -  (1/  (EXP  (x)  -  1))  -  (1  /  (EXP(x)  -  1)  ~  2) 

END  FUNCTION 

FUNCTION  Gx  (x,  ro)  STATIC 

'g(x)  =  pi (beta) /pO (beta)  -  ro  =  0  Root  Solved  by  Newton's  Method 
Gx  =  1  -  (1  /  x)  +  (1  /  (EXP(x)  -  1))  -  ro 
END  FUNCTION 

SUB  makeP  (p(),  e)  STATIC 
SHARED  DifMeth$ 

'Modified  From  K.  Mathews'  Original  makeP  SUB  by  G.  Sjoden 

'Exception:  LD  Pade  exponential  inserted  in  LN  model  by  K.  Mathews 
/ 

CONST  imax  =  3 

DIM  pp(-l  TO  imax)  AS  DOUBLE 
DIM  ee  AS  DOUBLE 

IF  LBOUND(p)  <>  -1  OR  UBOUND(p)  <>  imax  Tv’'  .  .  ;0? 

SELECT  CASE  D,.fMeth$ 

CASE  “LD" 
ee  =  CDBL(e) 

'Forward  Iteration:  Use  Pade  qpprox  to  EXP(-ee) 

'  (Allows  LN  algorithm  tc  give  LD  results) 
pp(-l)  =  (6  -  ee  *  2)  /  (6  +  ee  *  (4  ^  , 

p(-l)  =  pp(-l) 
pp ( 0 )  =  (1  -  pp ( -1 ) )  /  ee 

p ( 0 )  =  pp(0) 

FOR  i  =  1  TO  imax 

PP  ( i )  =  (1  -  i  *  pp  ( i  -  1 )  )  /  ee 

P ( i )  =  PP ( i ) 

NEXT  i 

CASE  ELSE  'All  other  Spatial  Quadratures 

SELECT  CASE  (e) 

CASE  0 

P<-1)  =  1 

FOP  i  =  0  TO  imax 
p(i)  =  1  /  (l  +  1) 

NEXT  i 
CASE  IS  >  0 

ee  =  CDBL(e) 
pp(-l)  =  EXP(-ee) 

p  (  - 1 )  =  pp ( - 1 ) 

pp(0)  =  (1  -  pp(-l) )  /  ee 

p ( 0 )  =  pp ( 0 ) 

ife  =  FIX(e)  'Fix  truncates  e  (does  NOT  round)  to  integer  value 

only 

SELECT  CASE  (ife  -  1) 

CASE  IS  <=  imax  'Partial  Forward  Iteration 

FOR  i  =  1  TO  (ife  -  1) 
pp ( i )  =  (1  -  i  *  pp ( i  -  1 ) )  /  ee 
p  (  i )  =  pp ( i ) 
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PRINT 

INPUT  “Choice?  (  DD,  LD,  SC,  LC,  LN,  SA,  LA,  EC,  MB,  EL  ):  “,  DifMeth$ 
DifMeth$  =  UCASE$ ( LEFT$ ( LTRIM$ { Di f Meth$ ) ,  2)) 

LOOP  UNTIL  INSTR ( " ELMBDDLDSCLCLNSALAEC " ,  DifMeth$)  AND  LEN ( DifMeth$ ) 

2 

SELECT  CASE  DifMeth$ 

CASE  “MB’1 

TotDi fMeth$  =  "  DD  LD  LC  SA  LA  EC  “ 

CASE  "EL" 

TotDi fMeth$  =  *  DD  LD  SC  LC  LN  SA  LA  EC  “ 

CASE  ELSE 

IF  INSTR ( Tot Dif Me th$,  DifMeth$)  THEN 
ELSE 

TotDi fMeth$  =  TotDifMeth$  +  LTRIM$ (DifMeth$  +  “  “) 

END  IF 
END  SELECT 

ELSEIF  UCASE$ (Response$)  =  “CS“  THEN 
TotDi fMeth$  =  “  NONE" 

ELSEIF  UCASE$ ( Response $ )  =  “N"  AND  TotDifMeth$  <>  "  NONE"  THEN 
IF  INSTR (TotDifMeth$,  "EC")  THEN 
CALL  SetowitchECtoLC 
LOCATE  12,  1 
FOR  i  =  1  TO  12 

PRINT  TAB (79) ;  “  " 

NEXT  i 

LOCATE  12,  1 
END  IF 

IF  Problemf i le$  =  ""  THEN 

LOCATE  12,  1:  PRINT  "WARNING:  No  Problem  File  is  Designated” 

SLEEP  (3) 

GOTO  Variables 
END  IF 

IF  UCASES (Out  f ile$ )  =  “OFF*  THEN 

LOCATE  12,  1:  PRINT  “WARNING:  No  Output  File  is  Designated" 

PRINT 

DO 

INPUT  "Proceed  ?  (  Y/N  ):  “,  a$ 
a$  =  UCASE$ (LEFTS (LTRIM$ (a$)  ,  1)) 

LOOP  UNTIL  INSTR (“YN“,  a$)  AND  LEN(a$)  =  1 
LOCATE  12,  1 

FOR  i  =  1  TO  10 

PRINT  TAB (79) ;  “  “ 

NEXT  i 

LOCATE  12,  1 

IF  a$  =  “N"  THEN  GOTO  Variables 
END  IF 
GOTO  Bottom 
END  IF 

GOTO  Variables 
Bottom : 

'  Open  tiles  as  Necessary 
IF  NOT  again%  THEN 
OPEN  Problemf ile$  FOR  INPUT  AS  # 1 
'check  for  correct  file: 

INPUT  #1,  Ident $ 

LOCATE  9,  38:  PRINT  Ident $:  LOCATE  12,  1 
DO 

INPUT  "Correct  File?  (  Y/N/Quit  ):  ",  a$ 
a$  =  UCASES  (LEFT$  (LTRIM$  (a$)  ,  D) 

LOOP  UNTIL  INSTR ( "YNQ" ,  a$)  AND  LEN (a$ )  =  1 
SELECT  CASE  a$ 

CASE  “N“ 
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CLOSE  #1 
Problemfile$  = 

Ident$  =  * “ 

GOTO  Variables 
CASE  “Q" 

END 

CASE  " Y " 

CLS 

CASE  ELSE 
BEEP 

PRINT  “ERROR:  Unsupported  choice  in  OpenProblemFile . " 

STOP 
END  SELECT 
END  IF 

Assign  First  String  toggles  and  Corresponding  Diff  Method 
ilotogl  =  1:  ihitogl  =  3 

DifMeLhS  =  LTRIM$ (MID$ (TotDifMeth$,  ilotogl,  ihitogl)) 

DifMeth$  =  UCASE$ (LEFT$ (LTRIM$ (DifMeth$) ,  2)) 

END  SUB 

SUB  BeginSets 

SHARED  Problemf ile$,  DifMeth$,  QuadMeth$,  fixup%,  FlfixS,  Itermax% 
SHARED  Change,  MomTol,  nk,  TotDifMeth$,  Outfile$,  Qoutf$,  TKfile$ 
SHARED  Qtkf $ ,  Ident$,  again%,  NewProblem%,  ResetAllFluxes$ 

SHARED  MomOIntFlag$,  MaxMomOrerr ,  MomlIntFlag$,  MaxMomlrerr 
SHARED  SwitchEC% ,  SwECtoLCSetPt 

t 

'  This  provides  the  user  with  input  values 

'  for  an  initial  run. 

/ 

Problemf ile$  =  “a : testl . IPT" 

DifMeth$  =  "DD" 

TotDifMeth$  =  "  NONE“ 

QuadMeth$  =  “S” 
fixup%  =  False 
Flf ix$  =  "NO" 

Itermax%  =  150 
Change  =  .00001# 

MomTol  =  IE- 10 
nk  =  8 
Ident$  =  "  " 
again%  =  False 
NewProblem%  =  True 
SwitchEC%  =  False 
SwECtoLCSetPt  =  .0005 
Reset A1 1 Fluxes $  =  "YES" 

Outfile$  =  "off" 

Qout  f $  =  "n" 

TKf ile$  =  "off" 

Qtkf  $  =  "n" 

MomO IntFlag$  =  “  0th:  CONSERVED":  MaxMomOrerr  =  0 
Moml IntFlag$  =  "  1st:  CONSERVED":  MaxMomlrerr  =  0 
END  SUB 

DEFINT  A-H,  J,  L-M,  O-Z 
FUNCTION  Choose  (prompt $)  STATIC 
CONST  False  =  0 
CONST  True  =  NOT  False 
DO 

PRINT  prompt$; 

INPUT  ch$ 

ch$  =  UCASES (LEFTS (LTRIM$ (ch$)  ,  1)) 
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LOOP  UNTIL  INSTR ( " YN10 " ,  ch$)  AND  LEN (ch$ )  =  1 
IF  ch$  =  " Y "  OR  ch$  =  "1"  THEN 
Choose  =  True 
ELSE 

Choose  =  False 
END  IF 

END  FUNCTION 
DEFSNG  H 

DEFDBL  A-G,  J,  L-M,  0-Z 
FUNCTION  DdxG  (x)  STATIC 

'  d/dx  g(x)  for  Newton's  Method  Root  Solving  beta  in  EC 
DdxG  =  1  /  x  /v  2  -  (1/  (EXP(x)  -  1))  -  (1  /  (EXP  (x)  -  1)  "  2) 

END  FUNCTION 

FUNCTION  Gx  (x,  ro)  STATIC 

'g(x)  =  pi (beta) /pO (beta)  -  ro  =  0  Root  Solved  by  Newton's  Method 
Gx  =  1  -  (1  /  x)  +  (1  /  (EXP (x)  -  1))  -  ro 
END  FUNCTION 

SUB  makeP  (p(),  e)  STATIC 
SHARED  DifMeth$ 

'Modified  From  K.  Mathews'  Original  makeP  SUB  by  G.  Sjoden 
'Exception:  LD  Pade  exponential  inserted  in  LN  model  by  K.  Mathews 

I 

CONST  imax  =  3 

DIM  pp(-l  TO  imax)  AS  DOUBLE 
DIM  ee  AS  DOUBLE 

IF  LBOUND(p)  <>  -1  OR  UBOUND(p)  <>  imax  THEN  STOP 
SELECT  CASE  DifMeth$ 

CASE  “LD" 
ee  =  CDBL(e) 

'Forward  Iteration:  Use  Pade  qpprox  to  EXP(-ee) 

'  (Allows  LN  algorithm  to  give  Ld  results) 
pp(-l)  =  (6  -  ee  *  2)  /  (6+ee*  (4  +  ee) ) 
p ( -  1 )  =  ppv-1) 

pp ( 0 )  =  (1  -  pp ( - 1 ) )  /  ee 

p ( 0 )  =  pp ( 0 ) 

FOR  i  =  1  TO  imax 

pp  ( i )  =  (1  -  i  *  pp  ( i  -  D)  /  ee 
p ( i )  =  pp ( i ) 

NEXT  i 

CASE  ELSE  'All  other  Spatial  Quadratures 

SELECT  CASE  (e) 

CASE  0 

p  (  -  1  )  =  1 

FOR  l  =  0  TO  imax 
p  (  i  )  =1/  (  i  +  1 ) 

NEXT  l 
CASE  IS  >  0 

ee  =  CDBL(e) 

pp(-l)  =  EXP ( -ee ) 

p  '  -1 )  =  pp( -1) 

pp ( 0 )  =  (1  -  pp ( -1 ) )  /  ee 

p ( 0 )  =  pp ( 0 ) 

ife  =  FIX(e)  'Fix  truncates  e  (does  NOT  round)  ro  integer  value 

only 

SELECT  CASE  (ife  -  1) 

CASE  IS  <=  imax  'Partial  Forward  Iteration 

FOR  i  =  1  TO  (ife  -  1) 
pp  ( i )  =  (1  -  i  *  PP  ( i  -  D)  /  ee 
p  (  i )  =  pp ( i ) 
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NEXT  i 

b#  =  1 . 7 857 142857 14286D-02  'Remainder  Backward  Iteration 

'p26  Starting  Value,  p27(ee)=l/56  Double  Precision 

For  Single  Precision,  compute  pl6  using  pl7(e)=l/36 
FOR  i  =  26  TO  5  STEP  -1 
b#  =  ( 1  -  ee  *  b# )  /  i 
NEXT  i 

FOR  i  =  4  TO  (ife  +  1)  STEP  -1 
b#  =  (1  -  ee  *  b#)  /  i 
p(i  -  1)  =  b# 

NEXT  i 

CASE  IS  >  imax 

FOR  i  =  1  TO  imax  'Total  Forward  Iteration 

PP  ( i )  =  (1  -  i  *  pp ( i  -  1))  /  ee 
P ( i )  =  pp(i) 

NEXT  i 
END  SELECT 
CASE  ELSE 

PRINT  :  PRINT  "Negative  argument  in  SUB  MakeP" 

PRINT  "Execution  Halted.":  STOP 
END  SELECT 
END  SELECT 
END  SUB 

FUNCTION  MAX  (x,  y)  STATIC 
IF  x  >  y  THEN 
MAX  =  x 
ELSE 

MAX  =  y 
END  IF 

END  FUNCTION 

FUNCTION  MIN  (x,  y)  STATIC 
IF  x  <  y  THEN 
MIN  =  x 
ELSE 

MIN  =  y 
END  IF 

END  FUNCTION 

SUB  MomOStepCheck  (FI,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu)  STATIC 
SHARED  MaxMomOrerr,  MomTol,  MomOIntFlag$ 

'Verify  Oth  Moment  Balance  of  the  Boltzmann  transport  equation  for 
each  direction  mu 

sumML  =  (Fr  -  FI  +  eps  *  Fa) 
sumMR  =  Sa  *  DeltaX  /  ABS(mu) 

IF  (sumML  <>  0)  OR  (sumMR  <>  0)  THEN 
Momrerr  =  ABS ( sumML  -  sumMR)  *  2  /  (ABS (sumML)  +  ABS ( sumMR) ) 

ELSE 

Momrerr  =  0 
END  IF 

IF  Momrerr  <  2  THEN  MaxMomOrerr  =  MAX (MaxMomOrerr ,  Momrerr) 

END  SUB 

SUB  MomlStepCheck  (FI,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu)  STATIC 
SHARED  MaxMomlrerr,  MomTol,  MomlIntFlag$ 

'Verify  1st  Moment  Balance  of  the  Boltzmann  transport  equation  for 
each  direction  mu 

sumML  =3*  (FI-  2*  Fa  -t-  Fr)  +  (eps  *  Fx) 
sumMR  =  Sx  *  Delta);  /  ABS  (mu) 

IF  (sumML  <>  0)  OR  (sumMR  <>  0)  THEN 
Momrerr  =  ABS ( sumML  -  sumMR)  *  2  /  (ABS ( sumML)  +  ABS ( sumMR) ) 
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ELSE 

Momrerr  =  0 
END  IF 

IF  Momrerr  <  2  THEN  MaxMomlrerr  =  MAX (MaxMomlrerr ,  Momrerr) 

END  SUB 

FUNCTION  ParseNthWord$  (a$,  n)  STATIC 

' for  a$  containing  words  separated  by  spaces 
'returns  the  n'th  word,  with  word  count  starting  at  1 
i  =  n 

'strip  leading  spaces 
b$  =  LTRIM$(a$) 

' remove  n  leading  words 
DO  WHILE  i  >  1  AND  b$  <>  “  “ 

'strip  first  word,  up  to  a  space  character 
k  =  INSTR ( b$ ,  "  ’) 

IF  k  THEN 

b$  =  RIGHT$(b$,  LEN(b$ )  -  k  +  1) 

END  IF 

'strip  leading  spaces 
b$  =  LTRIM$(b$) 
i  =  i  -  1 
LOOP 

'remove  trailing  spaces  and  words 
k  =  INSTR (b$,  "  “) 

IF  k  >  0  THEN  b$  =  LEFT$(b$,  k  -  1) 

ParseNthWord$  =  b$ 

END  FUNCTION 

SUB  PrtTK  (labl$,  vrbl ( ) ,  istart,  n)  STATIC 

I 

'  ->  NOTE  DEVICE  #3  currently  used  for  tk  file  output  <- 

'  Because  any  variable  beginning  with  H  will  be  single  precision, 

'  any  data  stored  using  this  routine  as  is  will  be  stored  in  single 

'  precision.  To  get  full  double  precision,  change  hvalue  below  back 
'  just  vrbl ( i ) . 

I 

PRINT  #3,  labl$ ;  ", “ 

J  =  0 

FOR  i  =  istart  TO  n 
hvalue  =  vrbl ( i ) 

J  =  J  +  1 

'  if  using  this  sub  to  read  tk  generated  data,  use  j<5 
'  because  tk  writes  5  datapoints  per  line,  comma  delimited 

'  j<4  is  used  here  to  accommodate  double  precision  numbers 

'  (if  used)  that  could  be  scrolled  across  a  carriage  return 
IF  (J  <  4)  THEN 

PRINT  #3,  hvalue;  ”,  * ; 

ELSE 

PRINT  #3,  hvalue;  ",  ” 

J  =  0 
END  IF 
NEXT  i 
PRINT  #3, 

END  SUB 

SUB  SetSwitchECtoLC  STATIC 
SHARED  SwECtoLCSetPt 

PRINT  "Set  Exponential  Characteristic  (EC)  Default  to  Linear 
Characteristic  (LC):* 

PRINT  "  Typical  Values" 

PRINT  "rhox-Sx/(3  Sa)  beta=-(b  dx)  [  0 . 5 ( 1 -rhox ) ] =pl (beta ) /pO ( beta ) " 
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PRINT  “  -or+  2E-6  +or-  1.2E-5  0.500001  or  0.499999" 

PRINT  “  -or+  2.5E-4  tor-  5.0E-4  0.50004167  or  0.49995833" 

PRINT  “  -or+  0.32  tor-  2.1  0.66  or  0.34  “ 

PRINT  “(beta  <=  2)  insures  (ABS(Sx)  <  Sa)  and  No  LC  Fixup  is  used" 
PRINT 

PRINT  "Current  Default  EC  to  LC  Method  when  ABS(beta)  <=  "; 

PRINT  USING  "##  ;  SwECtoLCSetPt 

INPUT  "  Change  This  ?  (Y/N) :  ",  a$ 
a$  =  UCASE5 (LEFT$ (LTRIM$ (a$)  ,  1)) 

IF  a$  =  "Y"  THEN 
DO 

INPUT  "New  Magnitude  of  beta  <=  2  for  EC  Default  to  LC  <=  ", 
SwECtoLCSetPt 

LOOP  UNTIL  SwECtoLCSetPt  <=  2! 

PRINT  USING  "  ##  .  ####^^«  ;  SwECtoLCSetPt; 

PRINT  "  Confirmed.":  SLEEP  (2) 

END  IF 
END  SUB 

SUB  StepDD  (Jin,  Jout,  Sa,  dx,  Sigma,  mu.  Fa)  STATIC 
SHARED  f ixup% 
eps  =  Sigma  *  dx  /  ABS(mu) 
a  =  eps  /  2 

Jout  =  (Jin  *  (1  -  a)  +  SGN(mu)  *  Sa  *  dx)  /  (1  +  a) 

'balance  equation  combined  with  diamond  difference  assumption 
IF  f ixup%  AND  (Jout  /  mu  <  0 )  THEN 

Jout  =  0  'negative  flux  fixup 

Fa  =  (Sa  +  ABS(Jin)  /  dx)  /  Sigma  'conservation  by  balance 
equat ion 
ELSE 

Fa  =  (Jin  +  Jout)  /  (2  *  mu)  'auxiliary  equation  (diamond 

diff  ) 

END  IF 

FI  =  Jin  /  mu 
Ft  =  Jout  /  mu 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  dx,  mu) 

END  SUB 

SUB  StepEC  (Sigma,  DeltaX,  mu,  FI,  Sa ,  Sx,  Fr,  Fa,  Fx)  STATIC 
SHARED  SwECtoLCSetPt,  SwitchEC%,  pi 
'Exponential  Characteristic  Spatial  Quadrature 
DIM  pe ( - 1  TO  3),  pb(-l  TO  3),  pbee(-l  TO  3) 
eps  =  Sigma  *  DeltaX  /  ABS(mu) 

GOSUB  Root SolveBeta 
'Compute  Walters  functions 

'  Note  p ( 3 )  values  are  not  required  in  this  scheme 
CALL  makeP(pe(),  eps) 

CALL  makeP(pb(),  ABS(beta)) 

CALL  makeP (pbee ( ) ,  ABS((beta  -  eps))) 

'If  beta  is  <0,  use  transformations  of  Walters  p()  functions  to  avoid 
'  iterations  of  negative  arguments  in  MakeP  SUB 
SELECT  CASE  beta 
CASE  IS  >=  0 
CASE  ELSE 

Be  =  EXP (ABS (beta) ) 
pO  =  pb ( 0 ) :  pi  =  pb(l):  p2  =  pb(2) 
pb(0)  =  pO  *  Be 
pb ( 1 )  =  (pO  -  pi)  *  Be 
pb(2)  =  (pO  -  2  *  pi  +  p2)  *  Be 
END  SELECT 

'Note  set  pbee ( 0 ) = (  pO(beta-eps)  exp(-eps)  ),  etc 

'If  (beta  -  eps)  is  <0,  again  use  transformations  of  Walters  functions, 
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'  and  note  pbee(l),  pbee(2)  are  not  currently  needed! 

SELECT  CASE  (beta  -  eps) 

CASE  IS  >=  0 

pbee { 0 )  =  pbee(O)  *  pe(-l) 

'NOT  NEEDED:  pbee ( 1 )  =  pbee ( 1 )  *  pe(-l) 

'  pbee (2)  =  pbee (2)  *  pe(-l) 

CASE  ELSE 

Be  =  EXP (ABS (beta  -  eps)) 
pO  =  poee(O) 

pbee(O)  =  (pO  *  Be)  *  pe(-l) 

'NOT  NEEDED:  pi  =  pbee ( 1 ) :  p2  =  pbee (2) 

'  pbee ( 1 )  =  ( (pO  -  pi)  *  Be)  *  pe(-l) 

'  pbee(2)  =  ( (pO  -  2  *  pi  +  p2 )  *  Be)  *  pe(-l) 

END  SELECT 

'Contributions  from  Entering  Flux  (FI) 

Frl  =  FI  *  pe(-l) 

Fal  =  FI  *  pe ( 0 ) 

Fxl  =  3  *  FI  *  eps  *  (pe (2 )  -  pe(l)) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

GOSUB  ECQuadrature 
CASE  0 

Frs  =  0 

Fas  =  0 

Fxs  =  0 

END  SELECT 
Fr  =  Frl  +  Frs 

Fa  =  Fal  +  Fas 

Fx  =  Fxl  +  Fxs 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

CALL  MomlStepCheck ( FI ,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu) 

Endof EC : 

EXIT  SUB 

RootSolveBeta : 

SELECT  CASE  Sa 
CASE  IS  <>  0 

rhox  =  Sx  /'  (3  *  Sa) 
ro  =  . 5  *  (1  -  rhox) 

CASE  ELSE 
rhox  =  0 
ro  =  .5 
END  SELECT 

'Determine  first  guess 
SELECT  CASE  ro 
CASE  IS  <  0 ! 
ro  =  .001 
betaO  =  -1000 
CASE  IS  >  1! 
ro  =  .999 
betaO  =  1000 
CASE  IS  =  .5 
betaO  =  0“ 

CASE  IS  >=  .77 

betaO  =  1  /  (1  -  ro) 

CASE  IS  >  .23 

betaO  =  12  /  pi  *  TAN((ro  -  .5)  *  pi) 

CASE  IS  <=  .23 

betaO  =  -1  /  ro 
END  SELECT 

'Rootsolve  for  beta  =  -(b  dx)  using  Newton's  Method  if  betaO  is  above 
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'  switchpoint  to  LC  (where  beta-0  and  closely  approximates  LC  with 
Sx<<Sa ) 

SELECT  CASE  ABS(betaO) 

CASE  0 

beta  =  betaO 
CASE  IS  >  SwECtoLCSetPt 
betatol  =  .00001 
GOSUB  LoopNewton 
CASE  ELSE 

SwitchEC%  =  True 

CALL  StepLC (Sigma,  DeltaX,  mu,  FI,  Sa ,  Sx,  Fr,  Fa,  Fx) 

GOTO  EndofEC 
END  SELECT 
RETURN 

LoopNewton : 
icount  =  0 
DO 

icount  =  icount  f  1 

beta  =  betaO  -  (Gx(beta0,  ro)  /  DdxG(betaO)) 

'PRINT  ro;  “  “;  icount;  “  ■;  betaO;  “  *  ;  beta 
correct!  =  ((ABS(beta  -  betaO)  /  (ABS (betaO )) )  <  betatol) 

IF  NOT  correct!  THEN  betaO  =  beta 
LOOP  UNTIL  (correct!)  OR  (icount  >=  1001) 

IF  (icount  >=  1001)  THEN 

PRINT  :  PRINT  "Newton  Iteration  Loop  in  StepEC  Exceeded  1000 
iterat ions " 

PRINT  "Execution  Halted.":  PRINT 
STOP 
END  IF 
RETURN 

ECQuadrature : 

'Compute  coefficient  a  (from  S(x)  =  a  exp(b  x) )  and  other  multipliers 

a  =  Sa  /  pb(0) 

cR  =  a  *  DeltaX  /  ABS (mu) 

Ca  =  cR  /  eps 

Cxi  =  3  *  Ca  *  beta 

Cx2  =  3  *  Ca  /  eps 

'Flux  contributions  from  sources 

Frs  =  cR  *  pbee(0) 

Fas  =  Ca  *  ( pb  <  0 )  -  pbee(0)) 

Fxs  =  (Cxi  *  ( pb ( 2 )  -  pb ( 1 ) ) )  +  (Cx2  *  (2  *  pb ( 0 )  -  (eps  +2)  * 
pbee ( 0 ) ) ) 

RETURN 
END  SUB 

SUB  StepLA  (Sigmk,  DeltaX,  mu,  FI,  Sa,  Sx,  Fr,  Fa,  Fx)  STATIC 
'Constrained  Linear  Characteristic  Spatial  Quadrature 
DIM  plt(-l  TO  3),  pe ( - 1  TO  3),  pt(-l  TO  3) 
eps  =  Sigma  *  DeltaX  /  ABS (mu) 

'Contributions  from  Entering  Flux  (FI) 

CALL  makeP^pet),  eps) 

Fr 1  =  FI  *  pe(-l) 

Fa  1  =  FI  *  pe ( 0 ) 

Fxl  =  3  *  FI  *  ( pe ( 0 )  -  2  *  pe ( 1 ) ) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

GOSUB  Const rainedL inear SourceQuadrature 
CASE  0 

Fr  s  =  0 
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Fas  =  0 
Fxs  =  0 
END  SELECT 
Fr  =  Frl  +  Frs 
Fa  =  Fal  +  Fas 
Fx  =  Fxl  +  Fxs 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

CALL  MomlStepCheck (FI ,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu) 

EXIT  SUB 

Const rainedLinearSourceQuadrature : 
a  =  Sa  *  DeltaX  /  mu 
x  =  Sx  *  DeltaX  /  mu 
SELECT  CASE  Sx 

CASE  -Sa  TO  Sa  '“WedgeLin" 

Frs  =  (a-x)  *pe(0)  +  2  *  x  *  pe ( 1 ) 

Fas  =  (a  -  x)  *  pe ( 1 )  +  x  *  pe(2) 

Fxs  =  3  *  (a-x)  *  (pe ( 1 )  -  pe(2))  +  x  *  (3  *  pe(2)  -  2  *  pe(3)) 
CASE  -3  *  Sa  TO  -Sa  '"WedgeNeg" 

rho  =  ABS(Sx)  /  (3  *  Sa) 

'PRINT  "WedgeNeg",  rho 
tau  =  (1  -  rho)  *  1.5 
CALL  makeP(pt(),  eps  *  tau) 

CALL  makeP(plt(),  eps  *  (1  -  tau)) 

Frs  =  2  *  a  *  plt(-l)  *  ( pt ( 0 )  -  pt(l)) 

Fas  =  2  *  a  *  (1  -  tau)  *  plt(O)  *  (pt ( 0 }  -  pt ( 1 ) ) 

Fas  =  Fas  +  a  *  tau  *  (2  *  pt ( 1 )  -  pt ( 2 ) ) 

Sum  =  (1  -  tau)  *  (pt ( 0 )  -  pt(l))  *  (plt(O)  -  2  *  (1  -  tau)  * 

plt(l)  ) 

Sum  =  Sum  -  tau  *  (1  -  2  *  tau)  *  pt ( 1 ) 

Sum  =  6  *  Sum  +  3  *  tau  *  (1  -  4  *  tau)  *  pt(2)  +  2  *  tau  ~  2  * 

pt  ( 3 ) 

Fxs  =  a  *  Sum 

CASE  Sa  TO  3  *  Sa  ' "WedgePos 

rho  =  ABS(Sx)  /  (3  *  Sa) 

'PRINT  "WedgePos",  rho 
tau  =  (1  -  rho)  *  1.5 
CALL  makeP(pt(),  eps  *  tau) 

Frs  =  2  *  a  *  pt (1) 

Fas  =  a  *  tau  *  pt ( 2 ) 

Fxs  =  a  *  tau  *  (3  *  pt ( 2 )  -  2  *  tau  *  pt ( 3 ) ) 

END  SELECT 
RETURN 
END  SUB 

SUB  StepLC  (Sigma,  DeltaX,  mu,  FI,  Sa,  Sx,  Fr,  Fa,  Fx)  STATIC 
'Linear  Characteristic  Spatial  Quadrature 
DIM  plt(-l  TO  3),  pe ( - 1  TO  3),  pt(-l  TO  3) 
eps  =  Sigma  *  DeltaX  /  ABS(mu) 

'Contributions  from  Entering  Flux  (FI) 

CALL  makeP(pe(),  eps) 

Frl  =  FI  *  pe(-l) 

Fal  =  FI  *  pe ( 0 ) 

Fxl  =  3  *  FI  *  ( pe ( 0 )  -  2  *  pe ( 1 ) ) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

a  =  Sa  *  DeltaX  /  mu 

x  =  Sx  *  DeltaX  /  mu 

Frs  =  (a  -  x)  *  pe ( 0 )  +  2  *  x  *  pe ( 1 ) 

Fas  =  (a  -  x)  *  pe ( 1 )  +  x  *  pe(2) 

Fxs  =3  *  (a-x)  *  (pe(l)  -  pe ( 2 ) )  +x*  (3  *pe(2)  - 
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2  *  pe ( 3 ) ) 


CASE  0 

Frs  =  0 
Fas  =  0 
Fxs  =  0 
END  SELECT 
Fr  =  Frl  +  Frs 
Fa  =  Fal  +  Fas 
Fx  =  Fxl  +  Fxs 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

CALL  MomlStepCheck ( FI ,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu) 

END  SUB 

SUB  StepLN  (Sigma,  DeltaX,  mu,  FI,  Sa,  Sx,  Fr,  Fa,  Fx)  STATIC 
'Linear  Nodal  Spatial  Quadrature 
DIM  pe ( - 1  TO  3) 

eps  =  Sigma  *  DeltaX  /  ABS(mu) 

'Contributions  from  Entering  Flux  (FI) 

CALL  makeP(pe(),  eps) 

Frl  =  FI  *  pe(-l) 

Fal  =  FI  *  pe ( 0 ) 

Fxl  =  3  *  FI  *  (pe ( 0 )  -  2  *  pe ( 1 ) ) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

GOSUB  LNLinearSourceQuadrature 
CASE  0 

Frs  =  0 
Fas  =  0 
Fxs  =  0 
END  SELECT 
Fr  =  Frl  +  F_  s 
Fa  =  Fal  +  Fas 
Fx  =  Fx!  +  Fxs 

CALL  M" .nOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

CALL  MomlStepCheck (FI ,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu) 

EXIT  SUB 

LNLinearSourceQuadrature : 
a  =  Sa  *  DeltaX  /  mu 
x  =  Sx  *  DeltaX  /  mu 
'Always  use  method  "WedgeLin" 

'Note  --  fixup  by  scalar  flux  rotation  may  be  used  in  main  program 
Frs  =  (a  -  x)  *  pe(0)  +  2  *  x  *  pe ( 1 ) 

Fas  =  (a-x)  *  pe ( 1 )  +x*pe(2) 

Fxs  =  3  *  (a-x)  *  (pe(l)  -  pe(2))  +  x  *  (3  *  pe(2)  -  2  *  pe ( 3 ) ) 
RETURN 
END  SUB 

SUB  StepSA  (Sigma,  DeltaX,  mu,  FI,  Sa,  Sx,  Fr,  Fa,  Fx)  STATIC 
'Constrained  Step  Characteristic  Spatial  Quadrature 
DIM  pi r ( - 1  TO  3),  pe(-l  TO  3),  pr(-l  TO  3) 
eps  =  Sigma  *  DeltaX  /  ABS(mu) 

'Contributions  from  Entering  Flux  (FI) 

CALL  makeP(pe(),  eps) 

Frl  =  FI  *  pe(-l) 

Fal  =  FI  *  pe ( 0 ) 

Fxl  =  3  *  FI  *  (pe(0)  -  2  *  pe ( 1 ) ) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

GOSUB  ConstrainedStepSourceQuadrature 
CASE  0 
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Frs  =  0 
Fas  =  0 
Fxs  =  0 
END  SELECT 
Fr  =  Frl  +  Frs 
Fa  =  Fal  +  Fas 
Fx  =  Fxl  +  Fxs 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

CALL  MomlStepCheck ( FI ,  Fr,  Fa,  Fx,  eps,  Sx,  DeltaX,  mu) 

EXIT  SUB 

Const rainedSt epSourceQuadrature : 
a  =  Sa  *  DeltaX  /  mu 
rho  =  ABS(Sx)  /  (3  *  Sa) 

SELECT  CASE  Sx 

CASE  0  TO  3  *  Sa  '"StepPos" 

CALL  makeP(plr(),  eps  *  (1  -  rho)) 

Frs  =  a  *  plr(O) 

Fas  =  a  *  (1  -  rho)  *  plr(l) 

Fxs  =  3  *  a  *  (1  -  rho)  *  (plr(l)  -  (1  -  rho)  *  plr(2)) 

CASE  -3  *  Sa  TO  0  '“StepNeg" 

CALL  makeP(plr(),  eps  *  (1  -  rho)) 

CALL  makeP(pr(),  eps  *  rho) 

Frs  =  a  *  pr(-l)  *  plr(O) 

Fas  =  a  *  (rho  *  pr ( 0 )  *  plr(O)  +  (1  -  rho)  *  plr(l)) 

Sum  =  (1  -  rho)  *  ((1  -  2  *  rho)  *  plr(l)  -  (1  -  rho)  *  plr(2)) 
Sum  =  Sum  +  rho  *  (pr(0)  -  2  *  rho  *  pr ( 1 ) )  *  plr(O) 

Fxs  =  3  *  a  *  Sum 
END  SELECT 
RETURN 
END  SUB 

SUB  StepSC  (Sigma,  DeltaX,  mu,  FI,  Sa,  Sx,  Fr,  Fa,  Fx)  STATIC 
'Step  Characteristic  Spatial  Quadrature 
DIM  plr ( - 1  TO  3),  pe(-l  TO  3),  pr(-l  TO  3) 
eps  =  Sigma  *  DeltaX  /  ABS(mu) 

'Contributions  from  Entering  Flux  (FI) 

CALL  makeP(pe(),  eps) 

Frl  =  FI  *  pe(-l) 

Fal  =  FI  *  pe ( 0 ) 

'Contributions  from  sources 
SELECT  CASE  Sa 
CASE  IS  >  0 

a  =  Sa  *  DeltaX  /  ABS(mu) 

Frs  =  a  *  pe(0) 

Fas  =  a  *  pe ( 1 ) 

CASE  0 

Frs  =  0 

Fas  =  0 

END  SELECT 
Fr  =  Frl  +  Frs 

Fa  =  Fal  +  Fas 

Fx  =  0 

CALL  MomOStepCheck ( FI ,  Fr,  Fa,  eps,  Sa,  DeltaX,  mu) 

END  SUB 
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